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INTRODUCTION 


Optical  methods  for  signal  processing  have  long  been  touted  as  playing  an 
important  role  in  the  future:  they  will  enable  complex  operations  to  be  performed  on 
large  arrays  of  data  at  a  very  rapid  rate.  This  prediction  is  based  on  the  inherent 

capability  of  optical  systems  to  operate  on  two-dimensional  (2-D)  data  planes  and  on  the 
ability  of  spherical  lenses  to  perform  the  Fourier  transform.  However,  the  promise  of 

optical  methods  to  rapidly  perform  signal-processing  tasks  remains  unfulfilled,  with 

certain  notable  exceptions  (e.g.,  processing  synthetic  aperture  radar  and  stellar  speckle 
interferometry  data).  The-e  are  several  reasons  for  this,  but  the  most  salient  are  the 
limitations  of  available  2-D  input/output  devices  (spatial  light  modulators  and  detector 
arrays),  the  fact  that  the  optical  phase  of  the  processed  signal  cannot  be  directly 
detected,  and  the  sensitivity  of  coherent  optical  systems  to  mechanical  disturbances  and 
speckle  noise. 

In  contrast  to  the  situation  for  2-D  optical  hardware,  signal  processor  technology 
for  temporal  (1-D)  signals  is  quite  advanced  in  capability  and  flexibility,  and  thus 

presents  the  interesting  prospect  of  applying  these  1-D  devices  to  2-D  signal  processing 
if  a  suitable  dimensional  transformation  can  be  employed.  In  effect,  this  would  allow 
the  rapid  parallel  processing  capability  to  be  "traded  off”  for  more  precise,  flexible,  and 
noise-immune  1-D  serial  processing  in  a  hybrid  system.  Several  dimensional 
transformations  are  available  for  deriving  1-D  signals  from  2-D  data  and  reconstructing 
processed  2-D  outputs,  the  most  familiar  being  the  television  raster.  But  another 
algorithm,  the  Radon  transform,  has  some  very  nice  mathematical  properties  that  make  it 
an  excellent  candidate  for  application  to  signal  processing.  These  properties  were  derived 
by  an  Austrian  mathematician,  Johann  Radon,  early  in  this  century,  and  the  transform 
bearing  his  name  has  become  well-known  in  recent  years  as  the  mathematical  basis  for 
medical  computed  tomography.  In  the  Radon  transform.  1-D  signals  are  derived  from 
2-D  input  data  by  "projection,"  i.e.,  integration  along  sets  of  parallel  lines.  The  2-D 
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signal  can  be  regenerated  by  "smearing"  and  summing  appropriately  filtered  1-D 


projections  back  over  the  2-D  space.  The  mathematical  properties  of  the  transform 
enable  2-D  signal-processing  operations  based  on  the  Fourier  transform  and/or 
convolution  operations  to  be  performed  by  means  of  1-D  operations  on  the  projections. 
Such  operations  include:  generation  of  complex  Fourier  transform.  Hartley  transform. 
Wigner  distribution  function,  general  2-D  filtering  and  correlation,  bandwidth 
compression,  spectrum  analysis,  and  cepstrum  analysis. 


THEORETICAL  INVESTIGATION  OF  THE  RADON  TRANSFORM 
APPLIED  TO  SIGNAL  PROCESSING 


In  his  original  development  of  the  mathematical  theory  of  the  transform.  Johann 
Radon  proved  two  theorems  that  have  been  the  basis  for  application  of  the  Radon 
transform  to  signal  processing:  the  central-slice  (or  projection-slice)  theorem  and  the 
filter  theorem.  They  demonstrate  that  2-D  Fourier  transforms  and  convolutions  can  be 
performed  by  1-D  operations  on  the  projection  data.  To  illustrate  mathematically,  a 
projection  of  a  2-D  function  f(r)  is  commonly  defined  by  a  linear  space- variant  integral 
transformation: 


Rj[f(r)]  =  Xf(p  ,0) 


•  po  oo 

d'r  fir)  <5(p  -  rrt). 
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(1) 


where  R*  denotes  the  Radon  transform  operator.  As  is  customary,  we  denote  scalar 
variables  and  ve- tors  by  normal-face  and  bold-face  characters,  respectively  The 

projection  Xf  is  a  function  of  two  variables:  the  radial  spatial  dimension  p  and  the 
azimuth  angle  <t> .  However,  all  of  the  operations  we  consider  operate  on  p  alone,  and 
therefore  we  can  consider  the  projections  Xf  to  be  1-D  functions  of  p  parameterized  by 
the  azimuth  angle  <t>  The  central-slice  theorem  states  that  the  Fourier  transform  of  the 
2-D  function  Hr)  is  obtained  by  performing  1-D  Fourier  transforms  of  each  projection 
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and  displaying  the  outputs  at  the  proper  radial  azimuths: 

SF  i  [f(r)]|  p  =fl„  =  £F,  [Xr(p,0)]  =  Af  (t\$) .  (2) 

The  geometry  of  the  Radon  transform  and  the  central-slice  theorem  are  shown  in  Figure 
1  of  Ref.  1.  The  filter  theorem  demonstrates  that  the  1-D  convolution  of  the 
projections  of  two  functions  at  the  same  azimuth  is  identical  to  the  projection  of  the 
2-D  convolution,  i.e., 

R2[f(r)  **  g ( r) ]  =  Rj[f(r)]  *  R2[g(r)],  (3) 

where  *  and  **  denote  1-D  and  2-D  convolution  respectively.  It  is  easy  to  see  that 
the  same  result  holds  for  correlation  operations  as  well.  The  processed  2-D  function 
may  be  reconstructed  using  any  of  several  algorithms  to  perform  the  inverse  Radon 
transform.2 

Our  analysis  of  2-D  operations  susceptible  to  solution  in  Radon  space  has  primarily 
exploited  these  two  theorems.  We  have  investigated  those  useful  signal-processing 
operations  that  can  be  decomposed  into  a  sequence  of  Fourier  transforms,  convolutions, 
and  other  achievable  1-D  and  2-D  operations  such  as  addition,  pointwise  multiplication, 
and  taking  logarithms.  Such  operations  include  Fourier  analysis  (computation  of  both 
the  power  spectrum  and  complex  transform),  the  Hartley  transform,  image  filtering  and 
correlation,  bandwidth  compression,  generation  of  the  Woodward  ambiguity  function  and 
the  Wigner  distribution  function,  some  spectrum  estimation  algorithms  (periodograms. 
Blackman-Tukey  analysis,  and  Yule-Walker  autoregressive  models),  and  the  cepstrum 
Work  by  other  authors  has  established1  that  the  Radon  transform  can  be  useful  for 
pattern  recognition  through  calculation  of  image  moments  and  the  Hough  transform 


CONSTRUCTION  OF  A  PRACTICAL  SYSTEM  FOR  2-D  SPECTRAL 
ANALYSIS  AND  IMAGE  FILTERING 

The  hybrid  system  constructed  to  perform  signal  processing  in  Radon  space  consists 
of  an  optical  scanner  (to  generate  the  Radon  transform  data).  1-D  signal  processors,  and 
a  computer-controlled  CRT  display.  The  optical  Radon  transformer  uses  a  laser  source, 
a  Bragg-cell  scanner,  and  anamorphic  optics  to  project  a  line-of-!ight  onto  a  2-D 
reflective  or  transmissive  object.  By  collecting  the  light  reflected  or  transmitted  by  the 
object  onto  a  detector,  a  signal  proportional  to  the  line  integral  of  the  reflectance  or 
transmittance  is  generated.  The  line-of-light  is  scanned  parallel  to  itself  by  the  Bragg 
cell  to  produce  a  temporal  signal  proportional  to  the  line-integral  projection  for  one 
azimuth  angle.  After  one  projection  is  generated,  the  azimuth  angle  is  changed  by  an 
image-rotating  prism.  Thus,  the  Radon  projections  are  generated  as  a  sequence  of 
temporal  electronic  signals.  For  obvious  reasons,  the  optical  Radon  transformer  is  called 
a  flying- line  scanner,  and  is  shown  schematically  in  Figure  2  of  Ref.  1.  Though  we 
had  originally  planned  to  demonstrate  Radon  transformation  at  video  rates  (30  frames/s), 
we  are  limited  by  the  rotation  rate  of  the  stepper  motor  for  the  image  rotator  to  about  5 
frames/s.  This  is  by  no  means  a  fundamental  limit  for  signal  processing  in  Radon 
space--optical  systems  have  been  built  to  rotate  images  at  75  frames/s  with  excellent 
stability  and  image  quality. 

After  derivation  of  the  projection  data,  signal  processing  can  be  performed  by  1-D 
electronic  or  hybrid  devices.  For  the  demonstration  of  2-D  spectrum  analysis  and 
Fourier  transformation,  we  implemented  the  chirp  transform  algorithm  with  surface 
acoustic  wave  dispersive  filters  to  produce  the  1-D  transform  of  the  temporal  input  data 
within  30  ms.  The  time-bandwidth  product  of  the  Fourier  transformer  is  only  50.  bui 
again  this  is  by  no  means  a  fundamental  limitation.  Filtering  of  the  1-D  signals  was 
performed  by  applying  the  projection  signal  to  one  port  of  a  monolithic  SAW  convolver 
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A  fast  ECL  function  generator  was  constructed  to  store  the  filter  function  to  be  applied 
to  the  other  port  of  the  SAW  convolver. 

To  construct  the  2-D  Fourier  transform  signal,  the  1-D  processed  signal  was 
displayed  in  the  proper  polar  format  on  a  computer-controlled  CRT.  The  results 
obtained  with  the  system  are  available  in  Ref.  4.  which  is  available  in  the  Appendix. 

It  had  been  our  intention  to  design  and  construct  a  custom  SAW  filter  to  perform 
the  filtering  operation  for  image  reconstruction  from  projections.  However,  the 
capabilities  of  the  available  photolithographic  facilities  were  not  adequate  for  the  task, 
and  instead  we  utilized  the  SAW  convolver  for  the  filtering  operation.  The  ECL 
function  generator  was  built  to  store  the  filter  function.  Recognizable  reconstructions 
were  derived  of  input  scenes  at  approximately  5  frames/s.  but  were  not  of  useful  quality 
for  two  reasons.  The  signal-to-noise  ratio  of  the  output  from  the  SAW  convolver  was 
not  adequate,  and  the  original  image  rotator  used  to  perform  the  inverse  Radon 
transform  exhibited  too  much  runout.  The  results  obtained  are  to  be  published  shortly. 

PROOF-OF-PRINCIPLE  EXPERIMENTS  FOR 
OTHER  PROCESSING  OPERATIONS 

Both  computer  simulations  and  demonstrations  in  hardware  were  performed  for  a 
number  of  the  2-D  processing  operations  listed  above,  including  Fourier  spectrum 
analysis,  complex  Fourier  transformation,  the  Hartley  transform,  data  compression, 
generation  of  the  Wigner  distribution  function.  power  spectrum  estimation' 
(periodograms.  the  Blackman-Tukey  algorithm,  and  the  Yule- Walker  autoregressive 
model),  and  the  cepstrum.  Most  of  these  results  have  been  reported  either  in  the  open 
literature1-*-'-'-7-*  or  by  presentation  at  technical  meetings.  Papers  dealing  with  the 
remaining  operations  are  in  preparation. 
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FEASIBILITY  OF  USING  THE  RADON  TRANSFORM 


FOR  3-D  DATA  PROCESSING 

An  architecture  for  a  3-D  image  processor  was  developed  prior  to  the 
commencement  of  the  contract  period,  and  so  the  work  in  this  program  concentrated  on 
investigation  of  materials  for  rapid  storage  and  retrieval  of  the  data  arrays.  The 
proposed  technique  utilizes  wavelength-multiplexed  storage  in  alkali-halide  crystals.  A 
theoretical  examination  of  data-storage  mechanisms  in  the  crystals  was  made  to  describe 
the  conditions  for  a  linear  relationship  between  exposure  intensity  (or  exposure  time) 
and  hole  depth.  The  two  data-storage  mechanisms  are  photochemical  holeburning  (PI1B) 
and  nonphotochemical  holeburning  (NPHB).  It  was  discovered  that  PHB  materials  do 
exhibit  the  necessary  linear  relationship,  but  NPHB  materials  do  not.  The  results  were 
reported  in  Ref.  9. 
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Application  of  the  Radon  transform  to  optical  production  of  the 
Wigner  distribution  function 
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Abstract.  The  Wigner  distribution  function  (WDF),  a  simultaneous  coordinate 
and  frequency  representation  of  a  signal,  has  properties  useful  in  pattern 
recognition.  Because  the  WDF  is  computationally  demanding,  its  use  is  not 
usually  appropriate  in  digital  processing.  Optical  schemes  have  been  developed 
to  compute  the  WDF  for  one-dimensional  (1  -D)  signals,  often  using  acousto¬ 
optic  signal  transducers.  Some  recent  work  has  demonstrated  the  computation 
of  two-dimensional  (2-D)  slices  of  the  four-dimensional  (4-D)  WDF  of  a  2-D 
input  transparency.  In  this  latter  case,  the  required  2-D  Fourier  transformation 
is  performed  by  coherent  optics.  We  demonstrate  that  computation  of  the  WDF 
of  real  2-D  signals  is  susceptible  to  Radon  transform  solution.  The  2-D  opera¬ 
tion  is  reduced  to  a  series  of  1  -D  operations  on  the  line-integral  projections.  The 
required  projection  data  are  produced  optically,  and  the  Fourier  transformation 
is  performed  by  efficient  1  -D  processors  (surface  acoustic  wave  filters)  by 
means  of  the  chirp-transform  algorithm.  The  resultant  output  gives  1  -D  slices 
through  the  4-D  WDF  nearly  in  real  time,  and  the  computation  is  not  restricted 
to  coherently  illuminated  transparencies.  This  approach  may  be  useful  in  dis¬ 
tinguishing  patterns  with  known  texture  direction.  The  optical  setup  is  easily 
modified  to  produce  the  cross-Wigner  distribution  function,  a  special  case  of  the 
complex,  or  windowed,  spectrogram. 

Keywords:  optical  pattern  recognition  optical  data  processing :  Wigner  distribution  func¬ 
tion  Radon  transform  surface  acoustic  wave  signal  processing. 
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1.  INTRODUCTION 

The  Wigner  distnbution  function!  WDF)  was  introduced  in  1932  as  a 
phase  space  representation  in  quantum  mechanics.1  Because  it  de- 
scnbes  a  signal  simultaneously  in  Fourier  reciprocal  variables,  it  has 
potential  applications  in  the  recognition  of  nonstationary  pat¬ 
terns.1  >  The  WDF  of  a  l-D  input  function  of  f(x)  is  a  2-D  function 
and  is  commonly  defined  as 

-» 

Wf(x<).u)  =  /  f(«o+ 
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=  k-2u  ^(*0  +  y)r(xo  ~  y).  ■  "> 

where  the  operator  represents  a  l-D  Fourier  transformation  of 
coordinate  x’  to  frequency  u.  In  the  2-D  case,  the  WDF  is  four- 
dimensional  and  is  defined  as 


Wf<r0.u)  =  /  /  f(r„  +  y)f*/r0  -  y) 

— oo  — oo 

(2) 

X  e~:’""-,'dV  . 


where  r0  and  r'  are  2-D  coordinate  vectors  and  u  is  a  2-D  spatial 
frequency  vector.  If  W^Xq.u)  is  evaluated  at  zero  frequency  and  a 
change  of  variables  is  performed,  the  WDF  becomes  an  autocon¬ 
volution.  Thus,  the  WDF  may  be  interpreted  as  a  generalized  auto- 
convolution  at  nonzero  frequency.4 

Several  authors5-'  have  reviewed  the  properties  of  the  WDF. 
including  some  aspects  that  make  it  suitable  for  implementation  by 
optical  processing  Most  importantly,  the  WDF  of  any  real  or  com¬ 
plex  function  is  real  (though  not  always  positive),  since  it  is  the 
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Fourier  transform  of  f[XQ  +  (x,'2)]f*[x0  —  (x7  2)],  which  is  Hermi- 
tian  with  respect  to  x'.  In  addition,  the  region  of  support  of  Wf(x0,u) 
is  identical  to  that  of  f(x)  in  both  the  coordinate  and  frequency 
domains. 

In  computing  the  WDF,  the  major  bottleneck  is  the  Fourier 
transformation.  In  the  case  of  a  1-D  (2-D)  input  signal,  a  full  1-D 
(2-D)  Fourier  transform  must  be  performed  for  each  value  of  the  1-D 
(2-D)  coordinate.  Were  this  to  be  done  digitally  in  the  case  of  a  1-D 
discrete  signal  of  n  samples,  it  would  require  n  multiplications  to 
produce  the  product  function.  A  total  of  nlog;n  multiplications  is 
needed  to  compute  the  subsequent  fast  Fourier  transform,  giving  a 
total  of  n  -+■  nlog,n  multiplications  per  point.  This  sequence  must  be 
evaluated  at  each  sample  in  the  sequence  [corresponding  to  each 
value  of  x  in  Eq.  (1)],  giving  a  total  of  n[n  ■+■  nlog.n]  —  n-log,n 
multiplications  to  compute  a  1-D  discrete  WDF.  For  a  2-D  nXn 
array.similarreasoningdemonstratesthatatotalofn:[n-  +  n-log,n‘] 
~  2n4log,n  multiplications  is  required.  The  motivation  to  find  opti¬ 
cal  processing  algorithms  is  quite  apparent,  especially  in  the  applica¬ 
tion  of  feature  detection  or  recognition,  due  to  the  large  quantity  of 
output  data. 

Several  schemes  have  been  developed  to  generate  the  2-D  WDF 
of  1-D  signals. ’_l'  Recent  work  by  Bamlerand  Glunder10  has  demon¬ 
strated  computation  of  2-D  slices  of  the  4-D  WDF  of  a  real-valued 
2-D  input  transparency.  The  product  function  was  produced  opti¬ 
cally  by  an  autocollimating  telescope,  and  the  Fourier  transforma¬ 
tion  was  performed  by  a  lens.  By  scanning  over  the  coordinates  of  the 
input  transparency,  all  2-D  slices  of  the  complete  4-D  WDF  can  be 
found. 

Computation  of  Fourier  transforms  is  also  susceptible  to  solution 
by  the  Radon  transform. 1 1 " 14  Data  of  dimension  m.  where  m  >2.  are 
reduced  to  1-D  by  integration  over  m  —  1  dimensions.  A  1-D  Fou¬ 
rier  transform  of  the  projection  data  yields  one  line  through  the 
origin  of  the  m-D  Fourier  transform.  Varying  the  projection  angle 
allows  building  up  the  complete  Fourier  transform.  This  procedure  is 
easily  adapted  to  computation  of  the  WDF  and  offers  advantages  in 
certain  applications. 

2.  RADON  TRANSFORM 

The  Radon  transform  has  received  much  attention  in  the  scientific 
community  since  the  invention  of  x-ray  computed  tomography  (CT) 
in  the  1960s.  It  has  been  used  in  the  fields  of  astronomy,  geology,  and 
nuclear  magnetic  reasonance. 1 1  Recently,  it  has  been  adapted  to 
feature  extraction  in  optical  data  processing.15  In  1917  Johann 
Radon  published14  the  mathematics  of  the  transform,  in  which  he 
proved  that  a  2-D  mathematical  function  can  be  reconstructed  from 
the  complete  set  of  its  line-integral  projections  The  basic  mathemat¬ 
ical  analysis  of  the  Radon  transform  is  straightforward  and  has  been 
considered  by  several  authors,1112  so  we  shall  only  touch  briefly  on 
the  main  points  relevant  to  2-D  Fourier  analysis. 

The  1-D  line-integral  projection  A.(p.<p)  of  a  2-D  function  f(r) 
along  azimuth  direction  <t>  ( relative  to  the  x-axis)  is  defined  as 


Fig.  1 .  Geometry  of  the  Radon  transform.  (Top)  Derivation  of  one  protec¬ 
tion  Mp.d)  by  line-integral  projection.  Line  integrals  are  evaluated  along  the 
azimuth  direction  (d  T  (»r/2)l  to  yield  the  protection  along  azimuth  direction 
Id)  The  unit  vector  rt  defines  the  azimuth  (d).  (Bottom)  Central-slice  theo¬ 
rem:  the  1  -0  Fourier  transform  of  a  line-integral  projection  yields  one  line 
through  the  2-0  Fourier  transform  of  the  original  2-D  function. 


-  /  / 


d2rf(r)6(  p  —  r-ri) 


The  projection  X  may  be  regarded  as  a  I-D  function  of  p,  param¬ 
etrized  by  i i>.  The  1-D  delta  function  in  the  integrand  reduces  the  area 
integral  to  a  line  integral  along  a  tine  normal  to  n  and  at  a  distance  p 
from  the  origin  (Fig.  1).  The  set  [X(p,<j>)]  for  all  azimuth  angles  d> 
constitutes  the  Radon  transform  of  f( r).  As  will  be  demonstrated,  the 
WDF  of  a  2-D  function  f(x)  may  be  computed  by  performing  opera¬ 
tions  on  the  line-integral  projections  of  an  easily  derived  2-D  func¬ 
tion,  reducing  the  2-D  computation  of  the  Fourier  transforms  for 
each  value  of  the  coordinate  vector  t  to  a  senes  of  1-D  operations. 
This  can  be  seen  if  a  1-D  Fourier  transform  of  a  line-integral 
projection  is  performed: 
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where  lowercase  letters  denote  functions  in  coordinate  space  and 
uppercase  letters  denote  the  Fourier  transforms,  and  p  is  the  2-D 
spatial  frequency  vector.  This  result  shows  that  the  1-D  Fourier 
transform  of  the  line-integral  projection  A(p,<d)of  the  2-D  function 
ft  n  welds  one  line  through  the  ongin  of  the  2-D  Fourier  transform  of 
!i  ri  i  Fig.  I ).  This  is  the  central-slice  or  projection-slice  theorem.  The 
advantage  of  using  the  Radon  transform  approach  to  2-D  Fourier 
transformation  results  from  the  ability  to  do  the  Fourier  transforma¬ 
tion  in  one  dimension  once  the  projections  are  available.  There  are 
several  efficient  1-D  processors  available  to  perform  the  Fourier 
transformation,  including  acousto-optic  cells,  charge-coupled  device 
i  CCD)  transversal  filters,  and  surface  acoustic  wave  (SAW)  disper¬ 
sive  delay  lines.  The  system  constructed  uses  SAW  delay  lines  in  the 
chirp-transform  algorithm,  as  will  be  discussed  shortly. 

3.  FLYING  LINE  SCANNER 

To  use  the  Radon  transform  to  compute  Fourier  transforms,  it  is  first 
necessary  to  produce  the  line-integral  projection  of  the  2-D  function. 
This  is  easily  done  optically  using  a  device  we  call  a  flying  line  scanner 
(Fig.  2.),  which  projects  a  line  of  light  onto  the  input  transparency. 
The  azimuth  of  the  line  of  light  can  be  selected  by  an  image  rotator, 
e.g..  a  dove  prism.  The  light  transmitted  through  the  transparency  is 
proportional  to  the  line  integral  of  the  intensity  transmission  along 
that  line.  An  acousto-optic  scanner  allows  the  line  of  light  to  be  swept 
perpendicular  to  itself  [i.e.,  varying  p  in  A(p,d>),  Eq.  (3)].  The  light 
transmitted  is  collected  by  the  photomultiplier  tube  (PMT),  whose 
output  current  in  time  is  proportional  to  A(p.<p).  Rotation  of  the 
dove  prism  varies  the  angle  <t>  and  allows  the  entire  set  [  A(  p ,<p)]  to  be 
collected. 


4.  SAW  CHIRP  FOURIER  TRANSFORM 

The  SAW  filter  is  an  acoustoelectric  device  that  can  be  designed  to 
have  one  of  a  wide  variety  of  impulse  responses.  It  consists  of  a 
piezoelectric  crystal  substrate  upon  which  is  deposited  a  pair  of 
conductive  interdigital  transducers!  Fig.  3).  A  rf  signal  applied  to  one 
transducer  produces  a  rf  field  between  the  fingers  of  the  transducer. 
This  field  distorts  the  crystal  piezoelectrically,  and  these  displace¬ 
ments  travel  along  the  crystal  surface  at  the  sound  velocity.  When  the 
acoustic  wave  reaches  the  second  transducer,  an  electric  field  is 
piezoelectrically  induced  in  the  conductor.  The  resulting  electric 
signal  is  the  convolution  of  the  input  signal  and  the  filter’s  impulse 
response  By  appropriate  design  of  the  interdigital  transducers,  the 
desired  response  may  be  obtained.17 

To  perform  Fourier  transformation,  three  filters  with  linear  FM 
impuUe  responses  are  required  for  the  chirp-transform  algorithm. 
The  impulse  response  of  a  linear  FM  filter  is 

h(t)  =  ei(<"o*a<>t  =  jiiuot  esiat2  _  (5j 


where  is  the  frequency  at  t  =  0  and  a  is  the  “chirp  rate." 

If  we  ignore  the  constant  frequency  a^,  a  signal  fj(t)  applied  to  a 
filter  of  impulse  response  hit)  =  e‘lat‘ will  produce  an  output  signal 
<o<‘): 

f0it)  =  f,lt)  V°'; 


J  dr  f,(r)  e10*' “  r>:  .  (6) 

— oo 


where  *  denotes  convolution.  Expanding  the  exponential  factor,  we 


Fig.  2.  Flying  line  scanner.  Collimated  He-Ne  laser  light  is  focused  onto  the 
Bragg  call  by  a  cylindrical  lens.  The  zero-order  diffraction  is  blocked  by  the 
filter,  and  the  first-order  beam  paaaes  through  to  the  image  rotator  The 
relay  optics  images  the  line  of  light  onto  the  transparency  f ( T) .  Application 
of  a  linear  FM  signal  to  the  Bragg  cell  scans  the  line  of  light  across  the 
transparency.  The  transmitted  light  is  collected  by  the  photomultiplier  tube. 
For  a  particular  azimuth  angle  p  selected  by  the  image  rotator,  the  PMT 
output  signal  in  time  is  proportional  to  the  line-integral  projection  A(p.p). 


Fig.  3.  Layout  of  a  simple  surface  acoustic  wave  filter.  An  impulse  ap¬ 
plied  to  one  transducer  produces  a  traveling  acoustic  wave  on  the  surface 
of  the  piezoelectric  substrate.  The  frequency  of  the  wave  is  determined  by 
the  spacing  of  the  fingers  in  the  interdigital  transducer  snd  the  amplitude 
by  the  smount  of  finger  overlap.  The  acoustic  wave  is  sampled  by  the 
output  transducer.  The  overall  filter  impulse  response  is  the  convolution 
of  the  responses  of  the  two  transducers.  For  linear  chirp  filters,  the 
response  to  an  impulsive  input  is  a  signal  varying  linearly  in  frequency 
over  time. 


obtain 

X 

f„(D  =  e‘ar2  f  dr  |fj(r)  eior~J  e--,<*rt  (7) 
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line  of  light  produced  by  the  flying  lino  (connor  (Fig.  2)  patooi  through  tha 
input  transparency  f(x.y)  3  f(7).  Tho  light  coiloctod  by  tha  PMT  pro¬ 
duces  a  time  signal  proportional  to  tho  lino- integral  projection  X|t,^|.  This 
signal  is  multiplied  by  an  upchirp  [h(t|  =»  s+'<‘t2]  and  convolved  with  a 
downchirp  [h(t|  =  e~,*t*].  The  output  is  demodulated,  giving  a  signal 
proportional  to  the  magnitude  of  the  1  -0  Fourier  transform  of  A(t.4).  This 
is  displayed  on  a  CRT  and  integrated  on  the  output  plane  (photographic 
film |  to  give  the  2-0  Fourier  power  spectrum. 


Identifying  t  as  wv  a  produces  the  equation 


30 


The  Fourier  transform  is  thus  obtained  in  three  steps: 

(1)  fj(t)  is  multiplied  by  e-la‘J  (premultiplication). 

(2)  This  product  is  convolved  with  a  filter  of  impulse  response  eiat‘. 

(3)  The  resultant  is  multiplied  by  e-lat  (postmultiplication). 


If  only  the  modulus  is  required,  the  postmultiplication  can  be 
deleted.  Of  course,  in  actuality,  the  filters  have  finite  time  windows  of 
width  T.  which  affect  the  limits  on  the  integrals  in  Eqs.  (6)  through 
(8).  and  overall  have  the  effect  of  convolving  the  result  with  a 
sinc(t/T)  function.  In  practice,  the  premultiply  and  postmultiply 
chirps  are  produced  by  applying  an  impulse  input  at  the  appropriate 
time  to  SAW  filters  whose  impulse  response  is  the  appropriate  chirp. 

A  Fourier  transformer  with  this  algorithm  was  constructed  using 
dispersive  filters  from  Andersen  Labs  (models  DS- 120- 10-20-25 1  A 
and  -25 IB).  The  time  dispersion  of  both  models  is  20  tis,  and  the 
bandwidth  is  10  MHz.  The  chirp  slopes  of  the  two  models  are  of 
opposite  sign.  The  time-bandwidth  product  of  the  system  (and  hence 
the  number  of  resolvable  spots  in  the  transform)  is  only  50.  but  with 
more  sophisticated  filters  the  time-bandwidth  product  could  be 
boosted  to  2000  or  more,  if  required. 

A  2-D  Fourier  spectrum  analyzer  was  constructed  using  the  flying 
line  scanner  to  produce  the  projection  and  the  SAW  filters  to  take  the 
transform  (Fig.  4  ).  The  transformed  signal  is  demodulated  and 
applied  to  the  z-axis  of  a  CRT.  For  each  projection,  this  gives  one 
line  through  the  2-D  Fourier  transform.  For  each  new  azimuth,  a 
new  line  is  written  on  the  CRT  and  displayed  on  the  output  plane  at 
the  proper  orientation  by  the  image-rotating  dove  prism.  Results  of 
the  Fourier  analysts  of  a  test  pattern  are  shown  in  Fig.  5.  Taking  the 


Fig.  5.  S  pact  rum  analysis  by  moans  of  tha  Radon  transform,  (a)  Input 
function,  (b)  Output  obtained  from  apparatus  of  Fig.  4.  Tho  fundamontal 
spatial  frequencies  of  the  fine  gratings  and  three  orders  from  the  coarse 
grating  are  visible. 


line-integral  projection  data  requires  10  us.  and  the  transform  data 
are  read  out  less  than  20  ns  later,  so  it  is  feasible  to  perform  the  full 
2-D  spectrum  analysis  at  video  rates  if  the  image  rotation  rate  is  900 
rpm,  requiring  a  prism  rotation  rate  of  450  rpm. 

5.  RADON  IMPLEMENTATION  OF  THE  4-D  WDF 

To  compute  the  4-D  WDF  of  a  2-D  real  function  t(r).  it  is  necessary 
to  form  the  product  function  t{r0  (r  2)]  t[r0  -(r  2)]  —  m(r0,r) 
for  all  values  of  r’  and  then  Fourier  transform  over  r\  We  can  apply 
the  Radon  transform  to  this  computation  in  the  following  manner 
First,  we  take  line-integral  projections  of  the  product  function 


Mp.r0.<p>  -//  dVm(r0.r)<5(p  -  r'-ri )  i9) 


The  geometry  of  the  projection  is  shown  in  Fig.  b  Taking  the  1-D 
Fourier  transform  of  A(p,r0.d>)  yields  [by  the  central-slice  theorem. 
Eq.  (4)]  one  line  through  the  2-D  s/fce  of  the  WDF  evaluated  at  r0  By 
rotating  the  azimuth  <t>,  we  can  build  up  the  2-D  slice  in  exact  analogy 
to  the  2-D  spectrum  analyzer  By  sampling  over  the  two  coordinate 
dimensions,  the  complete  4-D  WDF  can  be  computed.  The  geometry 
for  the  Radon  transform  calculation  of  one  line  through  the  WDF  is 
shown  in  Fig.  7. 
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Fig.  6.  line-integral  projection! of  theWigner diitnbution  function.  Tho 
integration  ia  made  over  the  line  p  =*  r'  ■  rt. 


Fig.  7.  Orientation  of  oneoutputlineof  theWDF.  Consider  a  line-integral 
projection  of  the  product  function  at  an  angle  4  to  the  x-axie  at  a  point 
(x„.0).  The  line  of  the  WDF  so  obtained  is  oriented  in  4-0  WDF  output 
space  aa  shown,  where  the  y-coordinate  axia  haa  been  ignored. 


light  from  the  flying  line  scanner  paaaea  through  the  beam  splitter  onto 
the  transparency  centered  at  ~Q  -T  (7'/2).  The  light  tranamitted  ia  re- 
focuaed  onto  thejransparency  by  the  lena- mirror  system,  but  is  now 
centered  at  7„  +  (7'/2).  The  output  ia  reflected  by  the  beam  splitter  onto 
the  PMT.  The  PMT  output  is  Fourier  transformed  by  the  SAW  filter  as 
before  and  yields  one  line  through  the  WDF  of  t(7). 


mirror  tilting  may  be  done  quite  rapidly.  Results  are  shown  in  Fig.  9. 

This  method  offers  an  advantage  over  that  of  Bamler  and  Gliinder 
in  some  applications.  Since  the  Fourier  transformation  is  not  optical, 
coherent  illumination  is  not  required  if  an  appropriate  scanning 
technique  is  used. 

6.  COMPUTATION  OF  THE  CROSS- WIGNER 
DISTRIBUTION  FUNCTION  AND  ITS  RELATION  TO 
THE  SLIDING-WINDOW  SPECTRUM 

The  sliding-window  spectrum  of  a  function  f(r)  windowed  by  a 
function  g(r)  is  defined  as10: 

90  30 

Sfg(r'.u)  =  j  j  f(r -j'jg'fr — ^-je-',nu'r  d;r  .  (10) 

—30  —30 

From  Eq.  (2).  we  can  define  a  cross-Wigner  distribution  function 
(CWDF)  to  be 


To  produce  the  line-integral  projection  of  the  product  function, 
the  technique  used  by  Bamler  and  Gliinder 10  was  adapted  as  pictured 
in  Fig.  8.  The  input  transparency  is  placed  in  the  flying  line  scanner 
with  the  optic  axis  passing  through  the  point  r0  of  the  transparency 
The  transmitted  light  is  collected  by  a  lens,  focused  on  a  mirror,  and 
reimaged  by  the  lens  back  on  the  transparency  The  doubly  transmu¬ 
ted  light  is  reflected  out  of  the  system  by  a  beam  splitter,  collected, 
and  detected  by  the  PMT.  The  PMT  output  current  is  proportional 
to  the  line  integral  of  m(f0,r').  As  the  flying  line  is  scanned  across  the 
transparency,  the  temporal  signal  out  of  the  PMT  is  proportional  to 
the  integral  of  m(r0,r')  for  different  values  of  r'.  The  signal  is  Fourier 
transformed,  yielding  one  line  through  the  WDF.  Other  values  of  r0 
may  be  interrogated  either  by  moving  the  transparency  relative  to  the 
optic  axis  or  by  tilting  the  mirror.  Using  a  galvanometer  scanner. 


sc  sc 

Wfg(r,u)  =  j  j  f|r  ~  •yjg'fr  — ^-je-;:nu'r  d:r  (111 


By  changing  variables  in  Eq.  ( 1 1 )  to  q  =  f  2.  we  obtain 


Wfg(r.u)  =2  J  J  f(r  ■+-  q)g*(r  -  q)  e  dq:  (12) 

— 30  — 30 

Assuming  a  symmetric  window  function  f<t(r)  =  g(—  r)]  and  using 
Eq.  ( 10),  we  find 


K 


K 
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Fig.  9  1-0  slices  through  the  4-0  WDF  of  two  2-0  objacta.  In  aach  casa.  tha  uppar  traca  it  tha  signal  from  tha  PMT.  representing  A(p.e)  [Eq.  IS)].  Tha 
lowar  traca  it  tha  output  of  tha  chirp  Fouriar  transformar.  Sinca  thara  waa  no  poatmultipty  chirp,  tha  magnituda  of  tha  transform  modulates  tha  carrier 
fraquancy.  In  (a)  andlb).  tha  object  is  a  grating  of  25%  duty  cycle  (25%  opaqua.  75%  transparent)  in  a  circular  aperture.  Tha  envelope  of  tha  upper  trace  is 
dua  to  tha  line-integral  protection  of  tha  circular  aperture.  In  (al.  tha  grating  is  positioned  with  tha  optic  axis  centered  on  an  opaque  grating  line  (defining  ~0 
in  Eq.  (2)].  The  components  of  tha  product  function  exactly  "overlay.  ''  and  tha  WDF  at  this  coordinata  is  dominated  by  tha  fundamental  frequency  of  the 
grating.  In(bl.  tha  object  has  bean  shifted  (varying  ~9)  so  that  an  opaqua  grating  tine  of  one  shifted  function  'overlays"  the  transparent  region  in  the  other 
shifted  function.  Hence  the  WDF  is  dominated  by  a  frequency  twice  that  of  tha  fundamental  of  the  grating  In  (cl  and(d).  theobiect  is  a  Fresnel  zone  plate, 
and  the  coordinata  displacement  is  normal  to  tha  scanning  line.  S  (lifting  ona  zone  plate  relative  to  tha  other  results  in  a  linear  moire  whose  spatial  frequency 
increases  linearly  with  increased  shift. 


Wfg(r,u)  =2  y  /  f(q  +  r)g*(q  -  r)  e-2,n(:u)'R  d2q 


=  2Sf,(2r.2u) 


Thus,  by  computing  the  CWDF  of  a  function  using  a  symmetric 
window,  we  can  find  a  scaled  version  of  the  sliding-window  spec¬ 
trum  This  is  useful  in  some  pattern-recognition  applications  where 
the  local  frequency  spectrum  is  of  value.1* 

Evaluation  of  the  CWDF  is  also  possible  using  the  Radon  trans¬ 
form.  The  setup  is  shown  in  Fig.  10.  It  is  similar  to  the  system  for 
finding  the  WDF  except  that  the  reflecting  telescope  arrangement 
has  been  replaced  with  a  second  lens  and  transparency  to  supply  the 


window  function.  As  before,  one  line  through  the  spectrum  is  calcu¬ 
lated  at  a  time.  In  cases  of  directional  texture,  this  will  result  in  a 
reduced  throughput  of  insignificant  data.  Results  for  an  Air  Force 
three-bar  chart  are  shown  in  Fig.  1 1. 

7.  CONCLUSIONS 

We  have  demonstrated  a  hybrid  optical  analog  electronics  processor 
that  can  rapidly  compute  l-D  lines  through  the  Wigner  distribution 
function  and  cross-Wigner  distribution  function  of  real-valued  2-D 
inputs.  In  certain  pattern-recognition  applications,  such  as  recogni¬ 
tion  and  classification  of  scenes  with  directional  texture,  this  tech¬ 
nique  offers  advantages  over  digital  processors  in  speed  and  over 
other  optical  processors  in  output  configuration. 
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I.  Introduction  and  Definitions 

Tradition. illy,  the  motivation  for  processing  signals  by  optical  means  is  due  primarily  to 
two  factors.  The  first  is  the  ability  of  coherent  optical  systems  using  spherical  elements  to 
perform  the  Fourier  transform,  while  the  second  is  the  inherent  capability  of  optical 
systems  to  operate  on  two-dimensional  (2-0)  data  planes.  For  2-D  signals  (e.g.  images), 
optical  processing  is  of  obvious  utility,  but  even  if  the  signals  to  be  processed  are  one- 
dimensional  (1-D),  optical  techniques  may  allow  parallel  processing  of  several  channels. 
The  increased  system  throughput  thus  obtained  may  make  optical  processing  attractive 
relative  to  more  precise  (but  as  yet  slower)  digital  electronic  technologies. 

The  main  thrust  of  research  in  optical  signal  processing  has  been  directed  at  applying 
either  or  both  of  the  capabilities  of  rapid  Fourier  transformation  and  parallelism.  However, 
there  are  problems  restricting  the  utility  of  optical  processing  that  are  well-known  to  those 
working  in  the  field  and  which  diminish  its  attractiveness  relative  to  digital  electronic 
processing.  Primary  among  these  are  the  limitations  of  available  2-D  input/output  devices 
(spatial  light  modulators  and  detector  arrays),  and  (for  coherent  systems)  speckle  noise. 
These  limitations  are  responsible  for  restricting  ,the  use  of  optical  processing  to  a  few 
applications  in  which  they  are  not  significant  (e.g.  off-line  synthetic  aperture  radar 
processing).  In  marked  contrast  to  the  situation  for  2-D  hardware,  signal -processor 
technology  for  temporal  (1-D)  signals  is  quite  advanced  in  capability  and  flexibility,  and 
hence  it  may  be  profitable  to  apply  that  1-D  technology  to  2-D  operations,  if  possible.  In 
effect,  this  would  allow  a  trade-off  between  rapid  parallel  processing  and  precise  serial 
processing  in  a  hybrid  system.  Several  algorithms  are  available  to  derive  1-D  signals  from 
a  2-D  input  and  reconstruct  the  2-D  processed  signal.  A  familiar  example  of  such  an 
operation  is  the  television  raster,  which  creates  a  1-D  temporal  signal  from  2-D  imagery  by 
scanning  and  rederives  the  2-0  image  by  stacking  segments  of  the  temporal  signal  Rhodes, 


1981b).  The  raster  transduction  was  used  in  optical  signal  processing  by  Thomas  (196b)  to 
generate  a  2-D  array  from  a  long  1-D  temporal  signal  to  use  as  input  to  a  2-D  optical 
processor.  Several  other  dimensional  transduction  operations  were  considered  by  Bartelt 
and  Lohmann  (1981).  Qie  that  is  becoming  more  famiiiar  can  be  called  a  tomographic 
transformation,  where  a  1-0  data  set  is  derived  from  a  2-D  signal  by  integration  along  sets 
of  parallel  lines.  The  relation  between  these  t^vp  sets  of  data  has  some  nice  mathematical 
properties  that  make  the  transformation  potentially  very  useful  in  both  analog  and  digital 
signal  processing. 

I.A.  History  and  Development 

The  mathematical  basis  for  the  tomographic  transformation  was  derived  in  1917  by 
Johann  Radon,  an  Austrian  mathematician.  Radon  proved  that  the  complete  set  of  1-D 
projections  of  continuous  2-D  or  3-D  functions  with  compact  support  contain  ail  of  the 
information  in  the  original  function.  The  projections  are  derived  by  integration  of  the  2-D 
function  over  sets  of  parallel  lines,  or  by  integration  of  the  3-D  function  over  parallel 
planes.  The  derivation  of  the  1  -D  projections  from  the  function  is  the  forward  Radon 
transform.  Radon  also  derived  expressions  for  reconstruction  of  the  function  from  its 
projections — the  inverse  Radon  transform.  Generalization  of  the  theory  has  made  it 
applicable  to  functions  of  higher  dimensionality  (John,  1955).  Another  development  was 
made  by  Cormack  (1963,  1964),  who  formulated  the  mathematical  expansion  of  projections 
into  circular  harmonics,  i.e.  a  discrete  angular  Fourier  series  representation  of  the 
projection  data. 

Radon  was  primarily  interested  in  using  projections  to  find  solutions  of  Fbisson's 
differential  equation  in  electrostatics,  but  his  work  has  been  applied  to  a  myriad  of 
scientific  disciplines  since  the  1950s,  including  crystallography,  radio  astronomy,  geophysics, 
nuclear  magnetic  resonance,  radiative  scattering,  and  diagnostic  radiology.  This  explosion 


Df  interest  is  evident  by  the  number  of  publications  on  the  subject,  especially  in  the  last  15 
years  or  so.  For  a  good  discussion  of  applications  and  an  extensive  bibliography,  see  Deans 
(1983).  No  doubt  the  application  of  the  Radon  transform  most  familiar  to  the  lay  public  is 
in  diagnostic  radiology.  The  new  fields  of  x-ray  computed  tomography  (CT),  emission 
computed  tomography  (ECT),  and  magnetic  resonance  imaging  (MR1),  which  enable  imaging 
of  cross-sectional  slices  of  the  body  of  a  patient  from  sets  of  projection  data,  have 
receiv  ?d  much  attention  in  the  popular  press.  Indeed,  the  medical  application  of  Radon's 
theory  is  the  source  of  its  now  familiar  name;  "tomography'  is  derived  from  the  Greek  word 
for  slice.  Each  of  these  new  medical  wonders  owes  its  existence  to  Johann  Radon  and  the 
subsequent  researchers  who  generalized  and  applied  the  mathematical  theory. 

In  each  of  the  applications  listed  above.  Radon's  mathematical  theory  is  used  to  solve 
an  inverse  problem,  where  the  source  function  is  mathematically  reconstructed  from  the 
proiection  data.  Of  course,  the  complete  infinite  set  of  projections  is  never  collected, 
making  it  impossible  to  uniquely  reconstruct  the  source  function;  only  some  'best’  estimate 
may  be  found.  We  shall  not  overly  concern  ourselves  here  with  such  niceties,  as  they  are 
somewhat  removed  from  the  purpose  at  hand  and  have  been  considered  at  length  elsewhere 
Rowland,  1979)  (Barrett  and  Swindell,  1981)  .  Rather,  we  wish  to  ir,  estigate  the  use  of 
the  Radon  transform  as  a  dime'  ional  transducer  In  signal  processing.  The  discrete  nature 
of  the  data  set  will  still  be  of  some  concern  to  us,  mainly  due  to  nonuniform  sampling  of 
Cartesian  space  by  the  transformation,  but  our  main  purpose  is  the  identification  of  signal 
processing  ooerations  that  are  possible  and  profitable  to  perform  via  a  tomographic 
transformation.  For  some  of  these,  the  processed  1-D  data  alone  may  be  sufficient  for  the 
task  at  hand,  but  often  it  will  be  desirable  to  reconstruct  the  processed  2-D  signal  from  the 
processed  protections  and  so  some  consideration  will  be  given  to  optical  methods  of 
generating  the  inverse  Radon  transform. 


1.8.  Basic  Theory 


In  the  literature,  there  are  several  extensive  mathematical  developments  or  the  theory 
or  the  Radon  transform,  e.g.  Helgason  (1980),  Deans  (1983),  and  Barrett  (1984)  . 
Consequently,  we  shall  keep  our  discussion  brief  and  emphasize  applicability  ratner  than 
completeness  or  mathematical  rigor.  Also,  we  shall  generally  restrict  our  treatment  to  the 
2-D  problem,  with  occasional  remarks  about  application  to  3-D  when  warranted. 

I.B.1.  Forward  Radon  Transform,  Projections 

Given  a  2-D  function  f(r)  =  f(x,y)  (as  is  common,  we  shall  denote  vectors  by  boldface 
characters),  a  single  projection  along  an  azimuth  angle  $  can  be  derived  by  integration 
along  all  lines  at  azimuth  $  *  ir/2.  The  one-dimensional  function  thus  generated  has  as 
independent  variable  the  perpendicular  distance  of  the  integration  line  from  the  origin. 
This  distance  is  the  magnitude  of  the  vector  p,  where  p  =  (p,$)  in  polar  coordinates.  It  is 

also  useful  to  define  a  unit  vector  n  =  =  (1,$)  =  [cos  $,  sin  *  ]  (n.b.  square  brackets 

denote  Cartesian  coordinates  and  parentheses  denote  polar  coordinates).  Naturally,  for 
each  set  of  integration  lines  at  different  angles  relative  to  the  x-axis,  a  different 
projection  is  derived.  \  common  notation  for  a  projection  is  X(p,$),  implying  that  X  is  a 
2-D  function.  But  since  ail  operations  on  the  projection  will  act  on  the  spatial  coordinate 
p  alone,  we  can  consider  the  projection  to  be  a  1-D  function  parametrized  by  the  azimuth 
angle  $.  Depending  on  one's  mathematical  preference,  X(p,$)  can  be  defined  in  a  number 
of  equivalent  ways.  For  example,  we  can  consider  a  projection  to  be  obtained  by 
integration  over  lines  parallel  to  the  y'-axis  in  a  system  of  coordinates  [  x',y’  ]  rotated  at 
angle  »  relative  to  the  original  [x,y]  axes.  However,  there  are  distinct  advantages 
obtained  by  defining  a  protection  as  a  2-D  integral  transform  whose  kernel  is  a  1 -D  Dirac 
delta  function  which  selects  the  projection  azimuth,  as  shown  in  Figure  1.  Consider  a 


prelection  azimutn  ;  defined  bv  the  polar  unit  vector  n  =  ;1,j).  We  wish  to  determine  the 
value  of  the  promotion  coordinate  p  that  will  be  influenced  by  a  point  in  the  2-0  function 
located  at  r  =  ( ;  r , ,  9 )  =  [r  cos  3,  r  sin  a].  As  is  apparent  from  Figure  1,  r  must  be  located 
on  the  line  normal  to  n  at  a  perpendicular  distance  from  the  origin  defined  by 

p  =  r  cos  (8-*)  =  (r  cos  3  cos  $  ♦  r  sin  9  sin  *  )=  r  *  n .  (1) 

Hence,  multiplying  f(r)  by  6(p  -  r  •  n)  collapses  the  area  integral  to  a  set  af  line  integrals 
for  the  azimuth  defined  by  n,  giving 
r®  f® 

X(P,*)  =  |  I  d*r  f(r)  6(p  -  r  *  n).  ,2} 

1  -<»J  -«• 

The  transformation  has  mapped  the  Cartesian  coordinates  [x,y]  to  a  new  system  (p,$), 
which  is  called  Radon  space.  We  have  a  choice  about  the  limits  on  the  new  coordinates.  If 
we  consider  p  to  be  bipolar  (-•  <_p  <  »),  then  X(p,$)  =  X  (.  -p,  ^  -n  ) .  We  may  therefore  limit 
t  to  the  region  iv0  <  i  <  t).  If  we  require  p  to  be  positive,  then  $  runs  over  2t  raaians. 

The  former  choice  is  usually  preferred,  since  it  simplifies  the  mathematical  development.  X 
plot  of  the  Radon  transform  in  (p,#)  space  (Figure  2)  is  termed  a  sinogram,  since  a  point  in 
Cartesian  space  maps  to  a  sinusoid  in  Radon  space.  From  eq.  (2),  it  is  easy  to  see  that  the 
Radon  transform  is  linear  and  space-variant.  It  is  often  convenient  to  express  the 
projection  operation  in  operator  notation,  e.g.  R?(f(r)]  =  X(p,$),  where  the  subscript  , 
denotes  that  the  function  being  transformed  is  two-dimensional. 

The  projection  operation  described  by  eq.  (2)  can  be  easily  extended  to  functions  of 
higher  dimensionality  (Barrett,  1984).  For  example,  a  1-0  projection  of  a  3-D  function  can 
be  obtained  by  integration  over  parallel  2-0  planes.  Hence  the  1-0  Dirac  delta  function  m 
eq.  (2)  now  reduces  the  volume  integral  to  a  planar  integral.  The  transform  collapses  the 
3-0  function  f(x,y ,z)  to  a  set  of  1-0  projections  (e.g.  \(p,},3)  )  parametrized  by  the  two 
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angles  denning  the  jmt  normal  to  the  planes  or  integration. 

I.B.2.  Central-Slice  Theorem 

Now  that  the  forward  Radon  transform  has  been  defined,  .ve  need  to  investigate  its 
properties  that  may  be  useful  for  signal  processing.  Foremost  of  these  is  the  central-slice 
theorem,  which  relates  the  fourier  transform  of  a  2-D  function  to  the  1-D  Fourier 
transforms  of  its  projections.  The  theorem  arises  because  the  kernel  of  the  Radon 
transform  is  a  Dirac  delta  function  of  the  scalar  product  of  the  conjugate  variables  r  ana  p, 
as  the  kernel  of  the  Fourier  transform  is  a  function  of  the  scaiar  product  of  conjugate 
variables  r  and  p.  As  is  customary,  we  define  the  Fourier  transform  of  a  2-D  function  f(r) 


f  •  f  * 

,  [ f( r) ]  3  F(p)  >  d*r  f(r)  e'2*1**'', 


where  jr.  is  the  2-D  Fourier  transform  operator  from  coordinate  r  =  [x,y]  to  spatial 
frequency  p  =  ''  y ,n  1  •  In  this  natation,  functions  denoted  by  a  lower-case  character  are 
the  coordinate-space  representation  (e.g.  f(r)),  while  the  corresponding  frequency-space 
representation  is  signified  by  the  upper-case  character  (e.g.  F( p ) ) .  If  we  perform  the  1-0 
Fourier  transform  of  the  projection  defined  by  eq.  (2),  we  obtain 


3  t  [kip,*;.]  =  A  ( v  ,  ?  )  =  I  dp  X(p,*>  e'2lTipv. 


Substitution  of  eq.  ;2)  into  eq.  (4)  yields 


A (v,$>  -  :  dp 


d*r  f(r)  5(p  -  r  •  n )  *  e*2’^. 


Exchanging  the  order  of  integration,  we  oOtain 


I, 


So  the  1-0  Fourier  transform  of  a  Radon  projection  at  azimuth  angle  $  relative  to  the  x- 
axis  yields  one  line  through  the  origin  of  the  2-D  Fourier  transform  of  the  function  f(r). 
This  line  (central  slice)  in  Fourier  space  is  oriented  at  the  same  azimuth  angle  but 
relative  to  the  5 -axis  (Figure  1).  The  central-slice  theorem  can  be  represented  in  operator 
notation  by. 

Rf  ;u) 

It  is  important  to  note  that  the  2-D  frequency-space  representation  generated  via  the 
Radon-Fourier  transform  has  a  sampling  density  in  Cartesian  space  that  falls  off  as  v  "* 
(Figure  3).  This  sampling  nonuniformity  must  be  compensated  whenever  a  Cartesian-space 
representation  is  derived  from  a  Radon-space  representation,  e.g.  for  display  of  the  2-D 
Fourier  transform,  or  (as  will  be  shown)  when  reconstructing  the  2-D  source  function  via 
the  inverse  Radon  transform.  Also  note  that  the  duality  of  coordinate-  and  Fourier-space 
representations  ensures  that  a  dual  to  the  central-slice  theorem  exists.  That  is,  the  inverse 
Fourier  transform  of  a  projection  in  Fourier  space  is  a  central-slice  of  the  coordinate  space 
representation  of  the  2-D  function. 

A  theorem  similar  in  nature  to  central -slice  relates  parallel  projections  weighted  by  a 
phase  factor  to  parallel,  rather  than  meridional,  lines  of  the  2-D  Fourier  transform  ^Farhat 


* 
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=  [Af  X  Agl  =  $  »"*  [Af  1  *  /,*“  [Agl 

=*f*Xg,  (13) 

where  the  common  coordinate  variables  have  been  suppressed.  In  words,  this  shows  that 
the  projection  of  a  2-D  convolution  of  two  functions  is  the  1-D  convolution  of  the 
projections  of  the  functions.  From  this  conclusion,  it  is  just  a  very  short  conceptual  hop  to 
the  realization  that  the  same  relationship  holds  for  2-D  correlations.  Thus,  we  now  have 
the  mathematical  capability  of  deriving  the  projection  of  a  2-D  filtering  or  correlation 
operation  simply  by  performing  1-D  filtering  or  correlation  of  the  projections  of  the  original 
functions.  This  is  a  very  powerful  result  and  holds  much  promise  for  application  to  optical 
processing. 

I.B.4.  Inverse  Radon  Transform 

Since  most  of  the  research  into  the  Radon  transform  has  been  directed  at  the  solution 
of  inverse  problems,  there  has  been  a  plethora  of  publications  devoted  to  the  inverse  Radon 
transform.  Therefore  we  shall  limit  our  mathematical  discussion  to  a  straightforward 
derivation  of  the  inverse  transform,  with  some  comments  made  about  algorithms  appropriate 
to  optical  reconstruction  methods.  Readers  interested  in  an  in-depth  mathematical 
development  should  consult  some  of  the  other  literature,  notably  Rowland  (1979),  Deans 
(1983),  and  Barrett  (1984). 

The  inverse  Radon  transform  is  most  easily  derived  by  applying  the  central -slice 
theorem  to  the  poiar  form  of  the  inverse  2-D  Fourier  transform: 


l[F(P)l  =  f(r)  =  {  d80  !  do  o  F(p)  e*2lMp’r  . 


Invoking  the  central-slice  theorem  [eq.  (7)],  we  set  p  =  nv,  o  =  v,  90  =  ♦,  and  F(p)  = 


F(ftu)  =  A(v,»)  in  eq.  (14),  yielding: 
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f(r)  =  j  d$  J  dv  |v|  A(v,$)  e*2iri\»n*r 
t,  J  — 

n 

<4  [  [M  A(v,*)]  l|p  _  .  (15) 

9 

This  is  one  form  of  the  inverse  Radon  transform.  In  words,  it  reconstructs  a  2-D  function 
f(r)  from  a  complete  set  of  projections  \(p,$)  by  the  following  steps: 

(1)  1-D  Fourier  transform  A(p,$),  yielding  A(v,$); 

(2)  multiply  by  |  v  |  ; 

(3)  inverse  1-D  Fourier  transform  the  product  [|v|  A(v ,♦ ) ]; 

(4)  smear  this  1-0  function  perpendicular  to  the  line  defined  by  p  =  r  *  n; 

(5)  sum  over  all  angles  4. 

Step  4  generates  a  2-D  function  from  the  1-D  projections  and  is  referred  to  as  ‘back- 
projection*  since  it  is  the  complementary  operation  to  projection.  Step  2  is  a  filtering 
operation  in  Fourier  space  to  correct  for  the  sampling  nonuniformity  of  the  transformation 
from  Cartesian  to  Radon  space  mentioned  previously. 

It  is  instructive  to  rearrange  the  steps  to  obtain  another  recipe  for  the  inverse 
transform.  Back-projection  and  sunmation  (steps  4  and  5)  may  be  performed  first  to 
generate  a  2-D  unfiltered  summation  image  (sometimes  called  a  ’layergram* ).  The  point 
spread  fuiction  of  the  layergram  has  been  shown  to  be  p(r)  =  | r | ” 1  (Peters,  1974),  wmch 
implies  a  transfer  function  [|r|_l  ]  =  jp|“l.  This  distortion  may  be  corrected  by 
filtering  in  2-D  with  transfer  function  |pj,  an  operation  commonly  known  as  ' rho-filtenng' 
(often,  albeit  imprecisely,  the  1-D  filter  |vj  in  step  2  is  also  referred  to  as  a  rho-filter). 

In  reality  of  course,  the  noise  dominant  at  high  spatial  frequencies  requires  either  filter  to 
be  rolled-off,  or  'apodized.'  Since  our  rationale  for  signal  processing  in  "<adon  space  was 
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to  avoid  unnecessary  2-D  operations,  .ve  shall  not  consider  implementations  of  the 
alternative  recipe.  Interested  readers  should  consult  Barrett  and  Swindell  ,1977,  1981)  or 
Barrett  (1984). 

We  can  also  express  the  inverse  Radon  transform  in  operator  notation  (Barrett,  1984), 
expanding  the  operator  Rt*‘  into  the  sequences: 

1  -  ^  1  1  M  ^ l 

=  jV‘  \o\fr  3,  ,  (16) 

where  Bt  is  the  operator  notation  for  back-projection. 

The  inverse  Radon  transform  algorithm  [eq.  (15)]  can  be  recast  into  a  more  concise 
form  by  invoking  the  filter  theorem  of  Fourier  transforms  to  create  a  convolution  of 
functions  instead  of  a  product  of  their  Fourier  transforms.  That  is, 

i  Ivj  A(  1  =  h(p)  *  x(p,?i,  .17) 

where  h(p)  =  j».-1  [  j  ]  is  the  filter  function  in  the  coordinate  space  representation. 

Lighthill  (1962)  showed  that  h(p)  a  _l  [  jvj  ]  a  ~2Tlpt,  where  the  singularity  at  the 

origin  requires  that  it  be  interpreted  as  a  generalized  function  which  has  a  Dirac  delta 
function  at  the  origin.  A  realizable  interpretation  is  (Gmitro  et  al.,  1980) 

h(p)  =  lim 
e  -0 

Note  that  h(p)  is  bipolar.  We  can  now  represent  the  inverse  Radon  transform  in  one 
equation,  with  the  important  proviso  that  the  true  nature  of  the  filter  function  be 
recognized : 


‘p7 


,P!  >  e 


Ip  <  e 


;i8) 


f(  r ) 


d*^(p,*)‘  pIJp  =  r.h- 


(19) 


ws 


IP#  7 


1*1 


r 


The  operations  required  to  implement  this  algorithm  for  the  inverse  rtadon  transform  are  the 
source  of  its  common  name,  filtered  back-projection.  Integration  of  the  convolution 
product  by  parts  yields  other  possible  expressions  for  filtered  back -project ion  (Barrett, 


1984): 


«'i  -  jp  f 


d*  [X"(p,4)  *  ln  bill  p  _  r  .  jj 


“2?  |  d*  [  ^  ' 
J  9 


ln  ipi]]p  =  r-S 


where  P  if  dx]  denotes  the  Cauchy  principal  value  of  the  integral,  and  the  primes  (e.g. 

)  represent  derivatives  of  the  function  with  respect  to  p.  Each  representation  of 
the  inverse  rtadon  transform  (eq.  (19-22)  j  requires  a  bipolar  filter  function,  a  fact  having 
important  consequences  for  optical  implementation.  Which  representation  is  optimum 
depends  strongly  on  the  limitations  of  the  signal  and  available  hardware.  For  instance,  the 
dynamic  range  of  the  1-0  filter  function  In  , p j  in  eqs.  (21-22)  is  much  less  than  that  of 
-p_t  or  P[p'1],  thus  reducing  the  dynamic  range  required  of  the  1-0  convolver  at  the  cost 
of  increased  noise  inherent  in  taking  the  second  derivative  of  the  proiection. 

An  alternative  development  of  the  reconstruction  problem  was  made  independently  by 
Cormack.  Though  not  as  straightforward  in  application  as  filtered  back-projection,  we  shall 
discuss  it  briefly  because  it  can  potentially  be  implemented  by  optical  methods  ^Ein-Gal, 
1974)  (Hansen  and  Coodman,  1978).  Cormack's  development  is  based  on  the  periodicity  in 
angle  of  every  physically  realizable  object,  i.e.  f(r,9}  =  f(r,d  *  2ir  i.  As  a  result,  f(r,9)  can 
be  expanded  in  a  discrete  Fourier  senes  of  angular  basis  functions,  which  are  called  circular 
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harmonics: 

as 

f(r,9)  =  Z  fn  (r)  e'n9  ,23) 

n=-® 

where 

fn  (f)  =  2V  [  d9  f(r,9)  e',n9.  (24) 

J  -at 

Cormack  expanded  the  projections  \(p,$)  in  the  same  manner  and  derived  the  space-variant 
transformation  between  these  two  representations.  The  transformation  can  be  made  space- 
invariant  via  a  Mellin  transform  (Casasent  and  Psaitis,  1977),  and  can  then  be  processed  by 
optical  methods  (Hofer,  1979)  (Hansen,  1981a,  1981b).  However,  the  Cormack 
reconstruction  algorithm  is  not  directly  applicable  to  our  task  at  hand,  so  we  shall  not 
consider  it  further. 

I.C.  Application  to  Optical  Signal  Processing 

To  summarize  the  mathematical  development,  we  have  demonstrated  that  the  classic 
2-0  signal -processing  operations  of  Fourier  transformation  and  convolution  (filtering)  can 
be  performed  via  the  equivalent  1-0  operations  on  the  rtadon  projections,  producing  central 
slices  of  the  2-0  Fourier  transform  or  projections  of  the  2-D  convolution.  Of  course,  there 
are  optical  methods  available  for  performing  these  2-D  operations  as  well.  Coherent 
computation  of  the  of  the  2-0  Fourier  transform  has  always  been  the  basis  of  optical  signal 
processing,  but  limitations  of  speckle  noise  and  performance  of  available  spatial  light 
modulacors  have  generally  restricted  application  to  static  film  transparencies  in  liquid 
gates.  3v  placing  the  input  in  the  front  focal  plane  of  the  transform  lens,  the  correct 
magnitude  and  phase  of  the  2-0  Fourier  transform  are  produced  in  the  back  focal  plane 
(limited  by  lens  aberrations).  However,  the  phase  of  the  transform  is  coded  in  the  relative 


pnases  of  cne  conerent  wavefront  at  the  various  locations  in  the  Fourier  plane.  Preserving 
this  phase  information  requires  a  very  precise  and  stable  optical  configuration,  and  square- 
law  detection  necessitates  heterodyne  techniques  to  decode  it.  Optical 
convolution/correlation  can  be  performed  by  spatial  filtering  in  the  Fourier  plane  or  by  a 
joint  transform  arrangement  (Weaver  and  Goodman,  1966)  (Rau,  1966).  Problems  still 
abound,  however.  The  stability  and  positioning  requirements  are  stricter  yet,  generation  of 
a  true  complex  (magnitude  and  phase)  spatial  filter  is  nontrivial,  and  deriving  the  phase  of  a 
complex  convolution  remains  difficult.  Incoherent  optics  avoids  the  specide  noise  problems, 
and  architectures  are  available  for  performing  Fourier  transformation  and  convolution 
(Rogers,  1977)  (Monahan  et  a!.,  1977),  but  representation  of  negative  quantities  requires  a 
bias  or  two  signal  channels. 

On  the  other  hand,  the  corresponding  1-0  operations  of  Fourier  transformation  and 
convolution  can  be  performed  readily  and  rapidly  by  devices  based  on  electronics,  acoustic 
interactions,  or  charge  transfer,  tiy  constructing  optical  systems  to  perform  the 
dimensional  transduction  to  and  from  Radon  space,  we  can  utilize  these  technologies  to 
'perform  the  corresponding  2-0  operation.  By  so  doing,  we  may  be  able  to  loosen  the 
constraints  on  signal  input  format  and  system  stability,  at  the  cost  of  some  processing 
parallelism.  The  resulting  hybrid  systems  can  emphasize  the  strengths  and  minimize  the 
weaknesses  of  each  technology.  If  the  optical  dimensional  transducers  and  the  1-0 
processors  are  fast  enough,  we  may  still  be  able  to  perform  the  complete  2-0  processing 
operation  at  a  usefully  rapid  rate,  e.g.  30  frames/second. 

I.C.1.  Optical  Radon  Transformer 

The  forward  Radon  transform  [eq.  (2)j  is  generated  by  integrating  the  input  function 
f(r)  along  the  set  of  lines  perpendicular  to  the  azimuth  $.  This  can  be  done  optically  in 
several  ways,  depending  on  the  format  of  the  input  data  and  the  type  of  signal  processor  to 


be  used.  Radon  projections  can  be  generated  as  temporal  data  by  scanning  the  input 
function  with  a  line  of  light  (usually  from  a  laser)  and  integrating  the  resultant  intensity  on 
a  detector.  This  method  is  suitable  for  transmissive  or  reflective  input  data.  At  one 
instant,  the  detector  signal  is  proportional  to  the  line  integral  of  input  transmittance  or 
reflectance.  Sweeping  the  line  of  light  perpendicular  to  itself  generates  a  temporal  signal 
proportional  to  one  line-integral  projection.  The  azimuth  of  the  scan  can  be  optically 
rotated  (e.g.  by  a  dove  prism)  to  sequentially  derive  the  complete  set  of  projection  data. 
For  obvious  reasons,  this  optical  Radon  transformer  is  called  a  flying-line  scanner,  and  is 
shown  schematically  in  Figure  4  (Easton  et  al.,  1984).  Since  the  light  transmitted  or 
reflected  by  the  2-0  input  is  integrated  on  the  detector,  speckle  noise  is  irrelevant,  and  a 
laser  can  be  usefully  employed  as  a  light  source.  Indeed,  the  coherence  of  the  laser 
becomes  an  advantage,  as  it  allows  the  use  of  a  fast  acousto-optic  beam  deflector,  or  a 
slower  and  cheaper  holographic  deflector  ('hologon  scanner').  The  technology  of  optical 
scanners  and  image  rotators  permits  a  system  to  be  built  capable  of  performing  Radon 
transforms  at  video  rates  with  video  resolution  (30  frames/sec,  500x500  points).  This 
would  require  scanning  500  azimuth  angles  with  500  resolvable  data  points  per  scan  every 
30  mS.  Acousto-optic  Bragg-cell  scanners  capable  of  resolving  more  than  1000  points  per 
10  uS  scan  have  been  reported  (Gottlieb  et  al.,  1983).  To  preserve  the  phase  of  the 
proiection,  the  temporal  center  of  the  flying  line  scan  must  intersect  the  image  rotation 
axis  each  time.  i.a.  the  optical  rotation  axis  of  the  prism  must  not  wobble.  Scanning  a  full 
projection  set  in  30  mS  requires  an  image  rotation  rate-of  180®/30  mS  =  900  KPM,  implying 
a  prism  rotation  rate  of  450  RPM.  Such  systems  have  been  constructed  and  demonstrated 
Cmitro  and  Gindi,  1985).  Indeed,  much  higher  rotation  rates  have  been  reported  .vnile 
preserving  holographic  image  quality  (Stetson  and  Elkins,  1977).  Radon  transformers  based 
on  the  flying-line  scanner  are  most  useful  for  2-0  signals  on  transparencies  (e.g.  movies/ 
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or  for  real  reflective  scenes. 

Projection  data  can  also  be  generated  by  'collapsing'  an  image  of  the  2-0  signal  onto 
a  linear  array  or  imaging  detector  with  anamorphic  optics  (Figure  5).  The  anamorphic 
optical  element  can  be  a  cylindrical  lens  (Gindi  and  Gmitro,  1984),  or  a  coherent  optical 
fiber  bundle  (Farhat  et  al.,  1983).  Alternatively,  if  an  N-element  1 -D  linear  array 
detector  can  be  obtained  with  an  aspect  ratio  of  N:1,  anamorphic  imaging  is  unnecessary. 
An  array  detector  samples  the  projection,  making  this  arrangement  especially  useful  if  the 
data  is  to  be  processed  digitally.  An  image  rotator  is  still  required  and  hence  the 
projections  are  again  generated  sequentially.  This  type  of  system  is  adaptable  to  naturally 
illuminated  scenes  or  to  self-luminous  signals,  as  from  a  CRT. 

I.C.2.  1-0  Signal  Processor  Technologies 

■Vs  was  demonstrated  in  eqs.  (19-22),  the  inverse  Radon  transformation  requires 
convolution  of  the  projection  data  with  a  bipolar  1 -D  filter  function.  Tlierefore  we  shall 
now  shift  gears  somewhat  to  investigate  the  types  and  capabilities  of  available  1-D  signal 
processors.  These  will  be  lumped  into  four  categories:  electronic  devices  (both  digital  and 
analog),  charge-transfer  devices  (mainly  CCDs),  acousto-electric  devices  .primarily  those 
based  on  surface  acoustic  waves,  or  SAWs),  and  acousto-optics  (AO).  In  the  first  case,  the 
Radon  transform  allows  direct  application  to  2-D  problems  of  the  very  technologies  that 
optical  methods  are  supposedly  competing  against  on  the  signal-processing  battlefield. 
I.C.2.a.  Electronic  Systems 

Electronic  systems  (analog  and  digital)  for  processing  temporal  signals  are  no  doubt 
familiar  to  the  reader.  They  can  be  as  simple  as  an  RC  filter  or  as  complex  as  a  digital 
supercomputer.  The  accuracy,  precision,  stability,  and  flexibility  of  electronics  are 
products  of  many  decades  of  theoretical  and  engineering  effort,  with  the  result  that 
electronic  systems  are  generally  preferred  for  signal-processing  applications.  This  is  the 
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target  it  .vhich  proponents  or  optical  signal  processing  must  aim,  but  it  is  moving  ahead  all 
the  time.  New  materials,  sucn  as  Ca.As,  and  new  fabrication  technologies,  such  as  x-ray 
litnography,  promise  further  improvements  in  packing  density,  speed,  and  cost  of  electronic 
devices.  Even  the  traditional  advantage  of  parallelism  offered  by  optical  processing  is 
fading,  as  new  algorithms  and  chip  architectures  are  adding  parallel  capability  to  the 
electronic  world. 

Electronic  signal  processing  is  generally  divided  into  analog  and  digital  domains,  each 
having  its  own  advantages  and  disadvantages.  Analog  processing  represents  signal 
amplitudes  by  proportional  voltages  that  can  be  added,  subtracted,  and  divided.  Some 
nonlinear  operations  (e.g.  thresholding)  are  easily  performed  as  well.  Analog  processing 
*ith  active  and  passive  components  can  be  fast,  with  bandwidths  reported  to  =  2  GHz  for 
silicon  uevicus  and  up  to  20  GHz  for  GaAs  (tlierman,  198S).  wore  complicated  operations 
e.g.  multiplication,  root  finding)  are  possible  with  special  analog  modules,  but  operation  is 
much  slower  and  suQject  to  severe  limitations  in  linearity,  stability,  and  precision.  For 
some  applications,  the  restrictions  can  be  eased  bv  using  the  analog  voltage  signal  to 
modulate  a  radio-frequency  (RF)  carrier.  RF  devices  capable  of  several  usefui  operations 
are  available,  including  multiplication,  phase  shifting,  and  phase  detection.  Though  still 
limited  in  linearity  and  stability,  these  devices  can  be  profitably  used  for  analog  signal 
processing. 

The  advantages  of  digital  systems  are  well  known- -probably  too  well  known  to  the 
optical  processing  community.  But  they  have  their  limitations  too,  lack  of  speed  and  large 
power  consumption  being  two  of  the  most  important.  Sampling  limits  system  bandwidth  and 
subjects  the  sampled  signal  to  aliasing.  A/D  and  D/A  conversions  may  have  to  trade  speed 
for  precision  and  dynamic  range.  Clock  rates  are  limited  to  *500  WHz  for  silicon-based 
logic.  However,  improvements  are  being  made  continuously.  For  instance,  the  increased 
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nobilitv  of  gallium  arsenide  charge  carriers  allows  clock  rates  up  to  several  GHz  with  lower 
power  consumption  ,Bierman(  1985).  Generally  trie  limited  disadvantages  of  digital 
processing  have  been  more  than  offset  by  its  inherent  noise  immunity  and  linearity.  An 
unlimited  variety  of  signal-processing  operations  are  amenable  to  solution  by  digital  means, 
and  new  special-purpose  hardware  promises  to  increase  speeds  dramatically.  The  very- 
high-speed  integrated  circuit  (VHSIC)  program  of  the  Department  of  Defense  is  stimulating 
the  design  and  production  of  new  devices,  such  as  the  Westinghouse  complex  arithmetic 
vector  processor,  which  can  pe-form  a  1024-point  16-bit  complex  Fourier  transfonn  in 
130  uS,  compute  one  point  of  a  256-element  16-bit  correlation  in  6  uS,  and  multiply  a  64x64 
16-bit  matrix  by  a  64-element  vector  in  35  uS  (Marr,  1982).  Digital  parallel  operation  is 
becoming  .more  economical  as  design  costs  drop  and  fabrication  yields  increase,  but  cost  is 
still  a  significant  limitation  for  such  devices  and  is  likely  to  remain  so. 

I.C.2.b.  Charge-Transfer  Devices 

Charge-transfer  devices  can  store  and  manipulate  packets  of  electronic  charge  using 
two  structurally  different  circuit  technologies.  The  older  'bucket-brigade'  device  is  a 
senes  of  MOS  transistors  and  capacitors,  where  the  charge  is  moved  between  capacitors  by 
alternate  switching  of  the  transistors.  These  have  been  largely  superseded  by  charge- 
coupled  devices  (CCDs),  where  minority  charge  carriers  are  stored  under  closely-spaced 
electrodes.  Charges  are  moved  to  detectors  at  the  edges  of  the  arrav  by  sequential  pulsing 
of  the  electrodes.  The  most  familiar  use  of  CCD  devices  has  been  as  1-D  and  2-0  optical 
detector  arrays,  where  the  amount  of  charge  in  a  detector  cell  is  proportional  to  the  photon 
flux.  Hr  '/ever,  it  is  also  possible  to  use  them  as  signal  processors,  where  the  sampled  data 
values  are  denoted  by  the  varying  amounts  of  charge  .  Bv  moving,  summing,  and  detecting 
the  charge  packets  >n  various  wavs,  a  variety  of  processing  operations  can  be  performed. 

The  resulting  devices  are  an  interesting  hybrid  of  ana'og  and  digital  qualities,  since  :ne 


amplitude  of  each  discrete  sample  is  a  continuous  variable.  CCDs  can  obviously  be  used  as 
delay  lines,  with  applications  to  signal  time  and  bandwidth  compression.  Tapped  delay  lines 
and  fixed-transversal  filters  can  be  constructed  by  spacing  nondestructive  charge  detectors 
along  the  charge  pathway  and  summing  the  tapped  signals  (Buss  et  al.,  1973)  (Beynon  and 
Lamb,  1980).  With  variable  weights,  the  filter  is  programmable.  Multiplying  adjacent 
tapped  signals  from  two  CCD  delay  lines  and  summing  the  products  allows  computation  of  a 
discrete  convolution.  The  useful  dynamic  range  of  these  CCD  devices  is  limited  by  the 
quantum  noise  floor  and  tne  saturation  level,  with  typical  specifications  of  60-70  dB  30  dB 
for  the  convolver).  The  bandwidth  of  the  CCD  devices  is  determined  by  the  analog 
electronics  and  the  sampling  clock  rate,  ranging  from  a  few  Hz  to  5  MHz. 

By  combining  the  CCD  devices  described  above,  a  wide  variety  of  1-D  signal- 
orocessing  operations  is  possible.  The  utility  of  fixed  and  programmable  CCD  transversal 
filters  and  of  the  CCD  convolver  for  signal  processing  is  obvious.  Using  two  or  three  filters 
«vit h  linear  FM  (or  chirp)  impulse  responses,  the  chirp  z-transform  algorithm  can  oe 
implemented  (Rabiner  et  al.,  1969).  This  algorithm  will  be  discussed  in  some  detail  later. 
CCD  spectrum  analyzers  using  the  chirp  z-transform  algorithm  have  been  demonstrated 
which  are  capable  of  computing  a  512-point  z-transform  at  a  5  MHz  sampling  rate. 

I.C.2.C.  Acousto-Electric  Devices 

Piezoelectric  materials  distort  when  placed  in  an  electric  field,  ana  also  they  generate 
a  field  when  mechanically  stressed.  By  applying  a  modulated  RF  electric  field  to  a 
piezoelectric  medium,  a  corresponding  acoustic  distortion  is  generated  which  can  be 
processed  and  detected.  This  acoustic  wave  propagates  in  the  medium  at  a  characteristic 
velocity  vs  *  io~3  c.  Thus,  the  acoustic  wavelengths  are  much  shorter  than  the 
electromagnetic  wavelengths,  allowing  signal  processing  devices  that  are  manv  wavelengths 
long  to  be  constructed  in  small  packages.  Components  based  on  acoustic  waves  in  bulk 
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materials,  sucn  is  the  quartz  oscillator  and  delay  line,  have  been  usee  *or  man-,  years. 
More  recently,  however,  much  attention  has  been  paid  to  using  acoustic  waves  on  the 
surface  or  a  medium  ^surface  acoustic  waves  or  SAWs;  due  to  their  accessibility.  Once  a 
wave  has  been  generated  on  the  surface  of  a  medium,  it  can  be  sampled  at  any  point  in  its 
journey  along  the  surface.  A  diagram  of  a  simple  SAW  device  is  shown  in  Figure  b.  -V  pair 
of  conductive  transducers  is  deposited  on  the  surface  of  the  piezoelectric  crystalline 
medium.  The  input  signal  (often  on  a  carrier)  is  applied  to  the  input  transducer,  consisting 
of  a  set  of  interleaved  ‘fingers'  connected  to  buss  bars.  The  field  distorts  the  medium 
piezoelectrically,  and  the  acoustic  wave  travels  along  the  surface  of  the  crystal  to  a  similar 
transducer  where  it  generates  an  electric  RF  signal. 

If  we  think  of  the  SAW  device  in  Figure  6  as  a  delay  line,  the  sampling  or  the  acoustic 
wave  by  the  output  transducer  is  a  tapping  and  summing  operation  performed  in  parallel  for 
many  points  in  the  acoustic  wave.  Hence,  the  Saw  device  is  another  example  of  a 
transversal  filter.  Variation  of  the  spacing  and  overlap  of  the  transducer  fingers  produces 
different  impulse  responses,  allowing  a  wide  variety  of  operations  to  be  performed.  The 
utility  of  SAW  filters  is  such  that  several  design  procedures  have  been  developed  Matthews, 
1977)  (Gerard,  1978),  and  the  filters  themselves  are  manufactured  bv  standard 
photolithographic  techniques  (Smith,  1978).  SAW  bandpass  filters  are  available  for  center 
carrier)  frequencies  from  10  MHz  to  2  GHz  and  bandwidths  from  <100  kHz  up  to  50^  of 
center  frequency  (Morgan,  1985).  The  noise-limited  dynamic  range  is  typically  70  dB, 
comparable  to  that  available  from  CCDs.  Indeed,  it  is  interesting  that  CCDs  and  SAW 
devices  are  so  complementary,  offering  similar  signal  processing  capability  over  a  wide 
range  of  input  frequencies  .Roberts,  1977). 

Linear  FM,  or  crnrp,  SAW  fiiters  are  easily  made  and  nave  found  wide  application  to 
radar  systems  Xlauder  et  ai . ,  19b0)  Gerard  et  a!.,  1973).  More  recently,  they  nave  been 
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employed  in  spectrum  analyzers  and  Fourier  transfo nners  (Jack  and  Paige,  19731  yjack  et 
ai 1980).  The  transducers  are  designed  such  that  the  impulse  response  of  the  filter  is  a 
signal  of  linearly  increasing  or  decreasing  frequency.  SAW  interdigital  chirp  filters  are 
limited  to  bandwidths  of  about  500  MHz,  dispersion  times  of  50  yS,  and  effective  time- 
bandwidth  products  of  about  1000  (Morgan,  1985).  Frequency  dispersion  can  also  be 
achieved  by  spacing  acoustic  reflectors  on  the  substrate.  These  so-called  reflective  array 
compressors  (RACs)  have  been  reported  with  bandwidths  to  180  MHz,  dispersion  times  to 
90  uS,  and  time-bandwidth  products  of  16,200  (Gerard,  et  al.,  1977). 

Other  useful  SAW  signal  processors  can  be  made  by  utilizing  the  nonlinear  response  of 
the  substrate  to  severe  distortions.  If  strong  acoustic  signals  are  applied  to  each  end  of  a 
substrate,  the  two  waves  will  interact  nonlineariy  to  generate  higher  harmonics.  The 
second  harmonic  of  the  carrier  frequency  contains  information  about  the  product  of  the  two 
signal  amplitudes.  Integration  of  the  second  harmonic  frequency  over  the  substrate  by  an 
area  electrode  produces  a  temporal  signal  proportional  to  the  convolution  of  the  input 
signals.  Since  second  harmonic  generation  is  inefficient,  the  convolution  signal  will  be 
weak,  typically  80  dB  below  the  input  signal  levels.  Even  so,  noise-limited  dynamic  ranges 
of  60  dB,  and  spurious-signal-limited  dynamic  ranges  of  30  dB  have  been  reported  (Ash, 
1978).  Acoustic  convolvers  are  available  commercially  with  time-bandwidth  products 
approaching  2000  (Morgan,  1985). 

I.C.2.d.  Acousto-optics 

Acousto-optic  processors  are  reviewed  in  detail  elsewhere  in  this  volume,  so  we  shall 
discuss  their  capabilities  only  briefly.  As  mentioned  above,  an  RF  electromagnetic  wave 
can  be  transformed  into  an  acoustic  wave  in  a  medium  via  the  piezoelectric  effect.  The 
variation  in  material  density  modulates  the  refractive  index,  producing  a  phase  grating 
which  can  diffract  light.  Devices  based  on  the  interaction  of  sound  and  light  have  long 


been  used  in  signal  processing  as  efficient  1-0  spatial  light  modulators  and  beam  deflectors 
vKorpel,  1981).  Developments  in  materials  and  architectures  in  the  'ast  15  years  or  so  have 
led  to  new  applications  for  bulk  A-0  devices  in  signal  processing,  including  time-integrating 
and  space-integrating  correlators/convolvers  (Berg  et  ai.,  1979)  (Rhodes,  1981aj 
(Abramovitz  et  at.,  1983),  Fourier  transformers  (Lee  et  ai.,  1982)  (Pancott  and  Reeve, 
1985),  and  generation  of  1-0  time- frequency  representations  (e.g.  the  Woodward  ambiguity 
function)  (Athale  et  al.,  1983)  (Casasent,  1983).  The  interaction  of  light  and  surface 
acoustic  waves  has  also  been  applied  to  various  signal  processing  operations  (Das  and  Ayub, 
1982)  (Casseday  et  al.,  1983).  Indeed,  AO  devices  and  SAW  devices  are  inherently 
compatible,  for  the  obvious  reason  that  the  processing  mechanism  is  so  similar.  Limits  on 
carrier  frequency,  bandwidth,  and  dispersion  time  are  comparable  for  both  types.  AO 
materials  support  carrier  frequencies  in  the  range  of  (1  MHz  <_va  <_1  GHz),  with 
bandwidths  of  up  to  500  MHz,  interaction  times  of  up  to  80  US,  and  time-bandwidth 
products  greater  than  10,000  (Berg  et  al.,  1979). 

I.C.3.  Optical  Implementation  of  the  Inverse  Radon  Transform 

Having  discussed  the  technologies  available  for  1-D  signal  processing,  we  are  now 
ready  to  describe  methods  for  reconstructing  the  2-D  processed  signal  from  the  1-D 
projections.  Two  mathematical  algorithms  for  reconstruction  have  already  been  discussed: 
filtered  back -projection  and  circular  harmonic  expansion.  As  stated,  the  latter  is  more 
complicated  to  implement  and  not  as  appropriate  for  signal  processing  applications,  and  so 
will  not  be  considered  further  here.  Interested  readers  should  consult  the  work  of  Hansen 
and  Goodman  (1978),  Hofer  (1979),  Hofer  and  Kupka  (1979),  and  Hansen  (1981a,  1981b). 

In  our  mathematical  development  of  filtered  back-projection,  we  stated  that  1-0 
filtering  can  be  performed  before  back-projection,  or  2-D  filtering  afterwards.  Optical 
reconstruction  systems  have  been  built  which  filter  in  2-D  (Peters,  1974),  but  3gain  we  are 
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more  concerned  with  application  or  1-0  technologies  to  the  problem.  Several  hvbrij. 
optical  systems  have  been  proposed  or  built  to  implement  1-D  filtered  back-protection,  and 
we  shall  give  a  brief  overview  of  those  systems  here,  Readers  desiring  more  detail  should 
consult  the  original  papers  or  the  review  articles  by  Barrett  and  Swindeil  (1977)  and  Gmitro 
et  al.  (1980).  To  lessen  problems  associated  with  coherent  noise,  these  systems  used 
incoherent  illumination.  However,  it  is  essential  to  recall  that  the  filtered  projection  is 
bipolar,  requiring  that  any  reconstruction  scheme  preserve  sign  information.  Because  of 
this  constraint,  systems  based  on  incoherent  optics  must  place  the  projection  signal  on  a 
bias  or  employ  two  signal  channels.  Neither  of  these  alternatives  is  desirable;  biased 
signals  reduce  the  contrast  of  the  reconstruction,  and  dual-channel  systems  are  subject  to 
differential  signal  errors. 

After  1-D  filtering,  the  algorithms  of  eqs.  (19-22)  require  two  more  steps:  back- 
projection  and  simmation.  Back-projection,  i.e.  generation  of  a  2-D  function  from  a  1-0 
projection  by  'smearing*  perpendicular  to  the  projection  azimuth,  has  been  demonstrated  by 
anamorphic  optics.  The  projection  is  written  on  the  face  of  a  1-D  display  device  (e.g.  a 
CRT  or  LED  array)  located  one  focal  length  from  a  cylindrical  lens,  and  imaged  onto  an 
integrating  2-D  detector  or  display  device.  As  this  operation  is  performed  for  each 
projection,  the  reconstructed  image  is  built-up  at  the  output  plane.  Any  integrating  2-D 
detector  can  be  used  for  summation  of  the  back-projections  (e.g.  photographic  film,  video 
camera,  or  human  eye  if  the  system  is  fast  enough). 

The  hybrid  optical-electronic  reconstruction  schemes  have  differed  greatly  in  detail 
and  degree  of  success.  The  system  of  Duinker  et  al.  (1978)  was  mostly  based  on  analog 
electronics,  with  only  filtering  performed  optically.  The  projections  were  displayed  in 
sequence  on  a  CRT  and  imaged  onto  two  area-weighted  optical  masks  representing  the 
positive  and  negative  parts  of  tne  filter  function.  The  images  of  the  protections  were 


swept  across  the  filter  masks  bv  electronic  deflection,  and  the  integrated  transmitted 
signals  electronically  subtracted  to  obtain  the  bipolar  temporal  filtered  signals,  back- 
projection  and  summation  were  performed  electronically.  Edholm  et  al.  (1978)  stored  the 
Radon  projections  on  film  in  sinogram  format  X(p,$).  A  filtered,  biased  sinogram  was 
generated  by  sandwiching  a  positive  image  of  X(p,$)  and  a  negative  image  of 
\(p,$)  *  n.(p),  where  h_(p)  represents  the  negative  part  of  the  filter  function  in  eq.  (lb). 
Back-projection  was  performed  for  each  line  of  the  sinogram  by  a  cylindrical  lens,  with 
summation  on  a  suitably  rotated  piece  of  photographic  film.  Despite  the  dynamic  range 
limitation  inherent  in  the  use  of  a  bias,  this  system  produced  some  good  reconstructions. 

Probably  the  most  successful  incoherent  optical  reconstruction  systems  synthesized  the 
required  filter  function  by  OTF  synthesis.  This  method  is  based  on  the  fact  that  the  OTF  is 
the  autocorrelation  of  the  pupil  function  (lohmann,  1977)  (Khodes,  1977)  (Rhodes  and 
Lohmann,  1978)  (Stoner,  1978).  Two  pupil  functions  are  calculated  for  which  the 
difference  of  the  autocorrelations  is  the  Fourier  transform  of  the  required  filter  point 
spread  function.  An  infinite  number  of  pairs  of  pupil  functions  are  theoretically  possible, 
with  the  optimum  choice  determined  by  system  requirements  such  as  light  throughput  or 
noise  considerations.  Since  the  required  positive  part  of  the  filter  psf  is  a  delta  function 
[eq.  (18)],  a  clear  pupil  in  the  positive  channel  is  appropriate.  Two  negative-channel 
pupils  successfully  demonstrated  are  the  so-called  Ronchi  pupil  (Barrett,  Creivenkamp  et 
al.,  1979),  and  a  logarithmic  phase  plate  (Barrett,  Chiu  et  al.,  1979).  The  envelope  of  the 
point  spread  function  of  either  pupil  falls  off  as  1/pJ,  as  required.  Optical  reconstruction 
systems  based  on  OTF  synthesis  include  the  drum  processor  (Gordon,  1977)  (Gmitro  et  al., 
1980),  the  loop  processor  (Greivenkamp  et  al.,  1981),  and  a  hybrid  digital -optical  system 
(Gmitro  et  al.,  1980).  An  example  of  image  reconstruction  with  the  loop  processor  is 


shown  in  Figure  7 


\  reconstruction  system  that  is  most  applicable  to  tomographic  signal  processing  tasks 
was  proposed  recently  by  Gmitro  and  Gindi  1,1985).  It  is  capable  of  performing  a  500  x  500 
point  reconstruction  of  projection  data  at  video  rates.  The  system,  depicted  in  Figure  8, 
implements  the  algorithm  of  eq.  (19).  Filtering  is  performed  by  a  space-integrating 
acousto-optic  convolver,  as  shown  in  Figure  9,  though  a  SAW  convolver  could  be  used  as 
suggested  in  section  I.C.2.C.  The  projection  data  are  stored  in  a  fast  digital  memory  and 
read  out  line-by-line  to  a  fast  O/A  converter.  The  analog  signal  modulates  an  KF  carrier 
and  is  then  impressed  on  a  Bragg  cell.  The  diffracted  light  is  Fourier  transformed  by  a  lens 
and  filtered  by  a  spatial  binary  transmission  mask.  The  diffracted  light  is  retransformed, 
collected  by  the  detector,  and  demodulated.  The  filtered  projection  is  displayed  on  a  CKT 
and  back-projected  by  a  cylindrical  lens.  Azimuth  selection  for  the  back-projection  is 
accomplished  by  an  image-rotating  prism,  and  the  2-D  image  is  collected  by  a  video  camera 
and  displayed  on  a  conventional  CKT.  The  image  data  are  read  out  rapidly  enough  for 
operation  at  video  rates  (30  reconstructed  frames/ second).  The  design  goal  is  to  process 
projections  at  video  rates  with  a  dynamic  range  of  12  bits,  implying  a  signal-to-noise  ratio 
of  about  4000.  Preliminary  results  are  presented  in  Figure  10. 
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II.  Applications 

II. A.  Operations  on  2-D  signals 

As  was  evident  from  our  mathematical  development,  the  application  of  the  Radon 
transform  to  signal  processing  primarily  exploits  the  central-slice  and  filter  theorems,  which 
allow  operations  based  on  Fourier  transforms  and/or  convolutions  to  be  performed  on  the 
1-D  projections.  Useful  operations  of  this  type  include  the  Fourier  transform  and  its 
relative,  the  Hartley  transform,  2-D  filtering,  some  pattern-recognition  algorithms, 
bandwidth  compression,  and  spectrum  estimation.  Some  of  these  operations  require  the 
flexibility  of  digital  operation  but  are  included  to  indicate  the  scope  of  application  of  Radon 
methods.  Since  application  of  projection  operations  to  signal  processing  is  a  field  that  has 
yet  to  be  fully  plowed,  much  of  our  treatment  will  deal  with  feasibility  rather  than  actual 
results. 

I  I.A.1 .  Fourier  Transformation 

Since  it  is  a  signal -processing  staple,  and  also  because  of  its  close  relationship  to  the 
Radon  transform  via  the  central-slice  theorem,  it  seems  natural  to  commence  our  discussion 
of  applications  with  2-D  Fourier  transformation.  After  having  been  generated  by  one  of 
the  systems  described  in  section  I.C.I.,  each  projection  is  Fourier  transformed  and  the 
result  is  displayed  in  the  polar  format  required  by  the  central -slice  theorem.  To  perform 
the  1-D  Fourier  transform,  we  introduce  the  chirp  transform  algorithm,  which  is  derived  by 
decomposing  the  Fourier  kernel: 

e-2irivt  _  x  e"iir(  St)1  x  e*'ir(j  *  3t)\  (25) 

The  three  complex  exponentials  are  quadratic  phase  terms  or  linear  FM  signals,  i.e.  the 
instantaneous  frequency  of  each  varies  linearly  with  time.  Such  signals  are  commonly 
called  chirps  by  the  radar  community.  The  factor  3,  with  dimensions  of  temooral 
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frequency,  nas  oeen  introduced  to  rationalize  the  units  or  the  exponent.  A  1-D  temporal 
Fourier  transform  can  now  be  written: 


Ft v )  =  dt  f(t)  e*2lTlvt 


=e*i’r(4)1  dt  [f(t)  j  x  e*'if(j  -  at)1 


=e-iir(  at)1 


x  [ [f(t)  e-<fT( at)2 ]  .  e-i-Tt at)* J  ^ 


Thus,  by  employing  three  temporal  chirp  signals  (one  with  positive  exponential  term,  or 
upchirp,  and  two  with  negative  terms,  or  downchirps),  the  Fourier  transform  of  f(t)  can  be 
implemented  in  three  steps: 

(1)  multiplication  of  f ( t)  by  a  downchirp; 

(2)  convolution  of  the  product  in  a  filter  with  an  upchirp  impulse  response; 


(3)  multiplication  by  a  downchirp. 

The  resulting  temporal  signal  is  a  scaled  version  of  the  Fourier  transform,  where  the 
frequency  is  related  to  the  output  temporal  coordinate  by  v  =  3*t.  The  pre-  and 
postmultiplication  chirp  signals  can  be  generated  by  applying  impulsive  inputs  to  filters  with 
upchirp  impulse  responses.  Note  that  this  analysis  has  assumed  that  the  chirp  signals  are 
complex  and  of  infinite  length.  If  only  the  power  spectrum  is  required,  the 
postmuitiplication  in  step  3  can  be  eliminated.  Because  of  the  order  of  operations,  this 
algorithm  is  usually  referred  to  as  the  M-C-M  chirp  transform,  for  multiply-convolve- 
multiply.  The  duality  of  multiplication  and  convolution  in  coordinate  and  Fourier  space 
imply  that  the  operations  can  be  exchanged  to  produce  a  second  arrangement  for  chirp 
transforms,  the  C-M-C  transform  (Jack  et  al.,  1980).  It  has  the  disadvantage  of  requiring 


three  filters  even  if  only  the  power  spectrum  is  required.  For  sampled  data,  Fourier 


■'  1ST. 
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transformation  is  equivalent  to  evaluation  of  the  z-transform  on  the  unit  circle.  The 
comparable  implementation  using  sampled  chirps  is  therefore  called  the  chirp  z-transform 
(Kabmer  et  al.,  1969). 

It  is  instructive  to  reconsider  coherent  optical  Fourier  transformation  in  light  of  the 
chirp  transform  algorithm.  Propagation  of  light  in  the  Fresnel  region  can  be  described  as 
convolution  of  the  wavefront  with  a  quadratic-phase  impulse  response,  and  the  action  of  a 
spherical  lens  on  a  wavefront  is  multiplication  by  a  quadratic  phase,  so  the  common  2*f 
coherent  Fourier  transformer  is  a  version  of  the  C-vt-C  chirp  algorithm.  An  optical  version 
of  M-C-M  is  also  possible  (Whitehouse,  1977). 

The  chirp  Fourier  or  z-transform  can  be  implemented  for  real  1-0  data  (as  would  be 
obtained  from  a  flying-line  scanner)  using  the  technologies  described  previously,  but  the 
analysis  differs  somewhat  from  that  given  in  eqs.  25-26.  A  basic  temporal  signal  filter  has 
a  real,  finite-length  impulse  response,  often  modulating  a  carrier.  For  example,  the  impulse 
response  of  a  SAW  chirp  filter  is  of  the  form 

h±(t)  =  A(t)  cos  <a,t  ±  ,  (27; 

where  A(t)  is  the  apodizing  function  of  the  filter,  a,  is  the  initial  carrier  frequency,  and  a 
is  the  'chirp  rate',  or  rate  of  change  of  the  instantaneous  frequency.  For  SAW  filters,  the 
carrier  frequency  a,  is  in  the  RF  region  (*15  -  300  MHz).  The  frequency  of  h*vt )  rises 
with  time,  so  this  function  is  again  called  an  upchirp.  Using  these  realizable  filters,  the 
chirp  Fourier  transform  may  still  be  implemented,  but  the  phase  of  the  transform  is  now 
determined  relative  to  the  phase  of  the  carrier  (Jac<  and  Paige,  1978).  The  recipe  for  the 
chirp  transform  becomes: 

(1)  premultiplication  by  a  downchirp; 

(2)  convolution  (filtering)  with  an  upchirp; 


1,3)  postmultiply  by  two  upcnirps  separately,  with  phase  difference  of  t/2; 

;4)  low-pass  filter  both  signals  from  step  3. 

The  complex  transform  is  thus  generated  as  two  parts  simultaneously.  The  signal  derived 
from  the  in-phase  chirp  of  step  3  is  the  real  part  of  the  complex  Fourier  transform,  or 
cosine  transform.  The  quadrature  signal  yields  the  sine  transform,  or  imaginary  part  of  the 
Fourier  transform.  Note  that  the  sign  of  the  slope  of  the  postmultiplication  chirp  differs 
for  the  realizable  algorithm  relative  to  that  for  complex  chirps.  This  is  due  to  double¬ 
sideband  multiplication  of  the  carrier-borne  signals,  which  yields  signals  at  the  sum  and 
difference  frequencies  of  the  carriers.  By  selecting  the  difference  frequency  sideband  with 
the  low-pass  filter,  the  operation  is  equivalent  to  postmultiplication  by  a  chirp  of  the 
opposite  sign.  The  output  temporal  signal  maps  linearly  to  frequency  with  constant  of 
proportionality  a.  Since  the  real  chirp  signals  are  apodized  by  A(t),  their  time-bandwidth 
product  (TBW)  is  finite,  thus  limiting  the  frequency  resolution  of  the  transformer.  The 
maximum  system  TBVV  is  one-fourth  the  TBW  of  the  convolution  chirp  (Ash,  1978).  It  should 
be  noted  that  the  SAW  chirp  transform  algoritfrn  can  also  be  implemented  for  complex  input 
data  by  premultiplying  the  imaginary  part  of  the  input  signal  by  a  chirp  in  quadrature  to  the 
real-part  premultiplication  chirp  (Jack  and  Paige,  1978).  Using  surface  acoustic  wave 
reflective  array  compressive  filters,  a  system  capable  of  transforming  signals  60  jS  long 
with  60  viHz  bandwidth  was  demonstrated  by  Gerard  et  al.  (1977).  SAW  chirp  Fourier 
transformers  are  faster  and  require  less  power  and  bulk  than  all-digital  systems,  but  are 
less  accurate. 

The  chirp  Fourier  transform  algorithm  can  be  implemented  with  40  devices  as  well. 

Hot z  1984)  and  Pancott  and  Reeve  (1985)  have  demonstrated  vf-C-M  transforms  using 
space-integrating  architectures  incorporating  two  Bragg  cells.  The  1-0  input  is  multiplied 
by  an  electronically-generated  chirp  signal  in  an  RF  mixer,  and  the  product  applied  to  one 


Bragg  ceil.  The  '1st  dirfraction  order  is  selected  and  imaged  on  the  second  Bragg  ceil, 
which  is  driven  by  the  same  electronic  chirp  signal.  The  -1st  diffraction  order  emerging 
from  the  second  ceil  is  selected,  integrated  on  a  detector,  and  demodulated.  Hotz  reports 
a  system  bandwidth  of  25  MHz  for  a  signal  duration  of  5  uS,  limited  by  the  capabilities  of 
the  AO  cells  and  by  problems  with  generating  the  proper  postmultiplication  chirp  slope. 
Such  a  system  has  similar  mechanical  stability  requirements  as  other  coherent  optical 
systems,  but  are  readily  applicable  to  signal  processing  in  Radon  space. 

M.A.I.a.  2-D  Power  Spectra 

Ticknor  et  ai.  (1984)  demonstrated  production  of  2-0  power  spectra  via  the  Radon 
transform  and  the  SAW  chirp  Fourier  transform.  Their  system  is  diagrammed  in  Figure  11. 
The  Radon  projection  of  a  2-0  transparency  f(r)  is  generated  by  a  Bragg-cell-driven 
flying-line  scanner.  One  projection  is  derived  tn  10  ■,»$.  Premultiplication  by  the  SAW  chirp 
is  performed  in  an  RF  mixer.  This  product  signal  is  applied  to  the  convolution  chirp  filter, 
whose  output  is  the  Fourier  transform  on  an  RF  carrier.  Since  the  phase  of  the  Fourier 
transform  is  not  required,  the  output  of  the  convolution  filter  is  detected  incoherently  with 
a  diode,  producing  a  unipolar  signal  proportional  to  the  squared-modulus  of  the  Fourier 
transform.  The  SAW  filters  used  had  time  dispersions  of  10  MHz  and  bandwidths  of  20  uS. 
Power  spectra  were  generated  by  the  system  within  28  uS  after  commencement  of  the 
flymg-line  scan.  The  spectra  were  20  uS  long  with  50  resolvable  frequencies.  By  the 
central-slice  theorem,  the  detected  signal  must  be  displayed  in  a  polar  format  to  generate 
one  line  through  the  2-0  power  spectrum.  However,  as  the  2-0  spectrum  is  built  up,  the 
polar  raster  oversamples  the  low  spatial  frequencies,  producing  a  displayed  time-averaged 
intensity  that  is  too  bright  in  the  center.  Mathematically,  this  problem  is  due  to  the 
sampling  nonumformity  of  the  Radon  transform,  and  is  corrected  by  rho-filtering,  i.e,  the 
central  slices  of  the  power  spectrum  are  multiplied  by  |v|  in  an  RF-mixer  before  detection. 


After  one  transform  slice  has  been  displayed,  the  prism  is  rotated  and  a  new  projection 
generated.  The  power  spectrum  of  that  projection  is  displayed  at  the  new  azimuth  on  the 
CRT.  Integration  of  the  result  can  be  done  on  film,  or  by  eye  if  the  system  is  fast  enough. 
System  speeds  up  to  5  frames/sec.  have  been  demonstrated,  limited  by  the  rotation  rate  of 
the  stepper  motor  driving  the  image  rotator  in  the  flying-line  scanner.  Results  for  a  2-D 
function  are  shown  in  Figure  12. 
li.A.I.b.  2-0  Complex  Fourier  Transforms 

The  same  group  (Easton  et  al.,  1985b)  added  a  post-multiplication  chirp  to  their  system 
to  generate  the  complex  Fourier  transform,  as  diagrammed  in  Figure  13.  The  time  delay  of 
the  post-multiplication  chirp  is  derived  from  a  digital  delay  generator  (1  ns  resolution).  To 
obtain  more  precise  time  delay,  the  phase  of  the  postmultiplication  chirp  can  be  varied  with 
a  continuously  adjustable  RF  phase  shifter.  The  postmultiplication  itself  occurs  in  an  RF 
phase  comparator,  which  generates  voltages  proportional  to  the  in-phase  and  quadrature 
products  of  two  input  signals.  The  in-phase  term  is  the  cosine  transform,  and  the 
quadrature  term  is  the  sine  transform.  Performance  of  the  complex  SAW  chirp  transformer 
is  shown  in  Figure  14. 

Rho-filtering  of  the  complex  transform  before  display  is  somewhat  more  difficult  than 
for  the  power  spectrum.  The  frequency  or  the  demodulated  signal  is  too  low  for 
multiplication  in  RF  mixers,  and  too  nigh  for  analog  multipliers.  An  integrated-circuit 
balanced  modulator  was  used  instead.  The  two  Dipolar  complex  Fourier  transform  signals 
were  then  biased  up  before  application  to  the  z-axis  of  the  CRT.  Results  are  shown  in 
Figure  15. 

Since  the  phase  of  the  transform  is  derived  from  the  time  differences  or  the  proiection 
signal  relative  to  the  chirps,  the  coherence  of  the  scanner  beam  is  immaterial.  This  method 
is  therefore  applicable  to  reflective  scenes  as  well  as  to  transparencies.  An  example  of 
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complex  Fourier  transformation  or  a  reflective  scene  is  shown  in  Figure  1b. 

Another  2-0  processing  situation  where  Raaon  space  Fourier  transformation  may  prove 
very  useful  is  with  spatial  light  modulators  whose  image  quality  is  relatively  poor. 
Recently,  there  has  been  much  interest  in  applying  an  inexpensive  liquid  crystal  television 
receiver  to  optical  processing  operations  (Liu  et  ai.,  1985)  (McEwan  et  al.,  1985).  The  poor 
phase  uniformity  of  the  LCTV  limits  its  utility  for  coherent  operations,  though  various  means 
have  been  suggested  for  improvement.  Again,  this  is  not  a  problem  when  used  as  input  for 
a  flying-line  scanner  (Easton  et  al.,  1985a).  Some  results  in  that  application  are  shown  in 
Figure  17. 

Farhat  et  al.  (1983)  also  demonstrated  2-D  complex  Fourier  transforms  via  Radon 
space  operations,  but  utilized  a  2-channei  incoherent  optical  correlator  to  generate  the 
1-0  transforms.  A  2-0  complex  signal  was  displayed  on  a  CRT  in  two  colors,  e.g.  real  part 
in  red,  imaginary  part  in  green.  The  image  was  rotated  bv  a  dove  prism,  spectrally  filtered 
to  separate  channels,  and  collapsed  to  1-0  by  two  coherent  optical  fiber  bundles.  The  real 
and  imaginary  1-0  signals  were  correlated  incoherently  with  a  fixed  cosine  and  sine 
reference  mask,  respectively.  The  1-0  correlator  outputs  represented  the  real  and 
imaginary  parts  of  the  1-0  Fourier  transform,  which  were  then  be  detected  and  displayed  in 
the  polar  raster.  The  system  is  fast,  but  also  suffers  from  the  familiar  limitations  on 
bandwidth  and  dynamic  range  common  to  other  geometrical-optics  incoherent  correlators 
Rogers,  1977). 

II.A.I.c.  Hartley  Transforms 

A  2-0  operation  that  is  receiving  some  attention  in  the  signal  processing  community  is 
the  Hartley  transform  (3racewell,  1983)  (Braceweil  et  al.,  1985).  For  a  2-0  function  f(ri, 
the  Hartley  transform  is  defined  as: 


where  cas(x)  3  cos(x)  +  sin^x).  The  kernel  or  the  Hartley  transform  Is  again  a  function  of 
the  scalar  product  of  conjugate  variables  and  is  in  fact  the  difference  of  the  real  and 
imaginary  parts  of  the  Fourier  kernel.  Being  purely  real,  the  Hartley  transform  may  be 
prererred  over  the  Fourier  transform  for  digital  computation,  since  the  storage  requirements 
could  be  halved.  Being  a  linear  combination  of  the  real  and  imaginary  parts  or  the  2-D 
Fourier  transform,  and  hence  of  the  1-D  Fourier  transform  or  the  Radon  projections,  the 
Hartley  transform  is  easily  implemented  in  Radon  space.  By  subtraction  of  the  real  and 
imaginary  outputs  of  the  SAW  chirp  transformer  with  a  simple  difference  amplifier,  the  1-D 
central  slices  of  the  Hartley  transform  are  genera:  d.  They  are  displayed  in  the  same 
fasmon  as  the  Fourier  transform. 

II.A.2.  Filtering  and  Correlations 

The  filter  theorem  demonstrates  that  a  projection  of  a  2-D  convolution  ^  correlation)  is 
the  convolution  (correlation)  of  the  corresponding  projections  of  the  2-D  functions.  Since 
devices  or  systems  exist  to  perform  1-D  convolutions  (SAW  devices,  CCD  convolvers,  and 
AO  convolution  systems),  it  is  feasible  to  perform  the  2-D  operations  in  Radon  space 
(Gmitro  et  al.,  1983).  With  a  fast  1-D  SAW  convolver,  such  an  operation  can  be  performed 
at  video  rates.  A  svstem  capaole  of  video-rate  2-D  convolution  or  filtering  is  depicted  in 
Figure  18.  The  projections  of  the  filter  function  may  be  generated  as  needed  from  a  2-D 
image  or  stored  in  digital  memory  and  read  out  through  a  fast  D/A  converter.  A  simulation 
of  2-D  nigh-pass  filtering  is  shown  in  Figure  19,  where  the  projections  were  generated 
optically,  the  1-D  convolutions  and  rho-riltenng  performed  in  a  digital  computer,  and  the 
back-projection  again  performed  optically. 


ir  the  prelections  or  the  filter  function  are  stored  in  an  addressaole  digital  memory,  as 
suggested  above,  we  have  the  capability  to  alter  the  impulse  response  of  the  2-0  filter  ov 
updating  its  digitally  stored  1-0  projections.  This  could  be  useful  if  filtering  a  noise  signal 
which  varies  over  time  and  would  enaole  the  application  of  1-D  adaptive  filtering  metnods 
to  2-0  situations.  For  example,  consider  a  signal  corrupted  by  noise.  An  adaptive  filter 
acts  on  noise  in  a  reference  channel  -.correlated  in  some  unknown  wav  with  the  noise  in  the 
signal  channel)  to  maximize  the  output  signal-to-noise  ratio.  This  is  accomplished  oy 
adjusting  the  filter's  impulse  response  to  minimize  an  appropriate  error  signal.  The  filter 
parameters  are  derived  from  correlations  between  the  signals  in  the  input  and  reference 
channels — operations  that  can  be  legitimately  performed  on  the  Radon  projections  of  2-0 
signals.  In  1-0,  the  technique  has  been  successfully  applied  to  a  number  of  problems,  e.g. 
telephone  echo  cancellation  (Gritton  and  lin,  1984),  electrocardiography,  and  antenna 
sidelobe  interference  jWidrow  et  al.,  1975).  To  the  knowledge  or  the  authors,  there  is  only 
one  demonstrated  example  of  2-0  adaptive  filtering.  Tao  and  iVeinhaus  ,1985)  applied 
adaptive  noise  cancellation  techniques  to  removed  of  periodic  signal-dependent  noise  in 
digital  imagery.  By  filtering  the  Radon  projections  with  1-0  updatable  stored  functions  in 
a  1-0  convolver,  these  adaptive  algorithms  can  be  implemented  while  avoiding  the 
limitations  of  available  2-D  hardware. 

II. A. 3.  Pattern  Recognition 

Some  very  useful  pattern  recognition  operations  can  be  proritaolv  performed  m  Radon 
space.  «Ve  have  already  demonstrated  generation  of  the  2-0  Fourier  power  spectrum. 
Cindi  and  Gmitro  ,1984)  have  used  optical  methods  to  rapidly  extract  integrated  features 
or  the  power  spectrum  from  the  Radon  projections.  They  have  also  demonstrated  the 
feasibility  of  evaluating  a  set  of  invariant  moments,  deriving  the  Hough  transform,  and 
rinding  the  convex  null  of  a  1-0  input  by  operations  -an  the  Radon  projections,  since  the 


first  three  operations  are  probably  of  most  interest,  we  shall  briefly  discuss  each. 

II. A. 3. a.  Fourier  Spectrum  Features 

Optical  computation  of  features  in  the  Fourier  power  spectrum  has  been  feasible  for 
some  years  and  has  been  applied  to  some  industrial  uses  (Casasent,  1981).  The  wedge-ring 
detector  was  developed  for  use  in  a  coherent  processor  to  compute  the  energy  in  the  power 
spectrum  in  discrete  segments  of  magnitude  and  orientation  of  spatial  frequency.  By 
manipulation  of  the  1-D  power  spectra  in  various  ways,  the  same  kind  of  Fourier  feature 
extraction  can  be  performed.  Integration  of  the  power  spectra  of  adjacent  projections 
produces  information  equivalent  to  that  from  the  wedge  segments.  Sampling  the  1-D 
spectra  and  integrating  over  projections  generates  information  from  discrete  spatial 
frequency  intervals,  corresponding  to  the  annular  segments  of  the  wedge- ring  detector. 
Results  from  a  computer  simulation  by  Gindi  and  Gmitro  (1984)  are  shown  in  Figure  20. 
II.A.3.b.  Image  Moments 

Two  decades  ago,  Hu  (1962)  described  a  system  of  linear  combinations  of  image 
moments  that  are  invariant  to  translation,  rotation,  and  scale  change.  Later,  Maitra  (1979) 
modified  the  system  to  include  invariance  to  relative  image  contrast.  Six  combinations  of 
ten  image  moments  m^  are  required,  where: 


m 


pq 


dx  dy  xP  yP  f(x,y)  . 


(29) 


The  ten  necessary  image  moments  are  m„,  mu,  nfij,,  mtl,  mJt,  mtl,  m,j,  m,„  m„,  and  m,,. 


Gindi  and  Gmitro  (1984)  demonstrated  that  the  ten  moments  can  be  computed  from  four 
projections  spaced  s/4  radians  apart.  The  ten  image  moments  and  the  linear  combinations 
can  be  rapidly  computed  by  digital  means  from  optically-generated  projections. 


1I.A.3.C.  Hough  Transform 

The  Hougn  transform  was  developed  as  a  technique  to  speed  detection  of  straignt  line 
segments  in  digital  imagery.  Edges  of  the  object  are  mapped  by  the  Hough  transform  to  a 
parameter  space,  wherein  peaks  indicate  the  presence  of  straight  lines  in  the  object. 
Deans  (1981)  described  the  close  similarity  between  the  Hough  and  Radon  transforms.  For 
binary  pictures,  in  fact,  they  are  identical.  Eichmann  and  Dong  (1983)  have  proposed  a 
coherent  system  to  generate  the  Hough  transform,  while  Gindi  and  Gmitro  (1984) 
demonstrated  that  1-D  filtering  of  Radon  projections  can  be  used  to  edge-enhance  a  2-D 
image  and  derive  the  Hough  transform  simultaneously.  Their  digital  simulations  of  the 
computation  of  the  Hough  transform  are  shown  in  Figure  20. 

I  I.A.4.  Image  Coding  and  Bandwidth  Compression 

The  potential  of  x-ray  tomography  in  medical  applications  led  to  investigation  of  the 
collected  data  required  to  obtain  good  image  quality  (Rowland,  1979).  In  turn,  this  has  led 
to  application  of  the  tomographic  transformation  to  reduce  image  storage  and  transmission 
requirements  while  maintaining  image  quality  (Mersereau  and  Oppenheim,  1974).  Since  only 
1-D  compression  operations  are  required  after  the  projections  are  collected,  rapid  coding  is 
possible.  To  date,  the  work  has  been  aimed  at  digital  compression  of  the  1-D  projections. 
Smith  and  Barrett  (1983)  truncated  and  quantized  the  Fourier  components  of  each 
projection  of  a  scene  to  reduce  the  data  from  8  bits/pixel  to  1.1  bits/pixel  while  retaining 
good  image  quality.  As  they  point  out,  the  approach  works  very  well  with  rectilinear 
scenes,  since  significant  Fourier  components  will  predominate  in  a  limited  number  of 
projections.  Fraser  et  al.  (1985)  investigated  the  effect  of  gross  reduction  of  the  number 
of  projections  used,  as  well  as  quantization  effects  of  various  spatial  frequency  ranges. 
Using  256x256  8-bit  images,  they  obtained  good  image  quality  with  as  few  as  0.86 


bits/ pixel. 


It. A. 5.  Spectrum  Estimation 

In  temporal  signal  processing,  the  estimation  or  frequencies  of  a  signal  buried  in 
uncorrelated  noise  is  a  classic  problem  (Robinson,  1982).  Averaging  and  modeling 
techniques  have  been  developed  appropriate  for  distinguishing  various  types  of  signals  (Kay 
and  Marple,  1981).  Most  are  based  on  Fourier  transform  and/or  correlation  operations  and 
are  hence  adaptable  to  operation  in  Radon  space  for  2-0  signals.  Traditional  methods 
incorporating  averaging  operations,  such  as  the  periodogram  and  the  Blackman-Tukey 
spectrum  estimate,  are  most  useful  for  detecting  the  presence  of  sinusoidal  signals.  To 
compute  the  periodogram,  windowed  segments  of  the  1-D  input  are  sampled  and  padded 
with  zeros.  The  size  of  the  data  window  determines  the  frequency  resolution  of  the 
periodogram.  The  power  spectra  of  the  segments  are  computed  and  averaged.  Since  the 
noise  is  uncorrelated,  the  signal  spectrum  should  dominate  in  the  periodogram.  This 
approach  has  become  popular  since  the  invention  of  the  FFT  algorithm.  2-0  periodograms 
are  used  in  a  similar  manner  for  spatial  signals  (Dudgeon  and  Mersereau,  1984).  For  2-D 
signals,  optical  processing  techniques  can  be  used  to  estimate  the  spectrum.  Indeed,  one  of 
the  success  stories  of  optical  processing,  Labeyrie  stellar  speckle  interferometry,  generates 
a  form  of  2-D  periodogram  where  the  signal  segmentation  is  over  time  rather  than  over 
space.  Computation  of  the  traditional  periodogram  is  readily  adaptable  to  Radon  space 
implementation.  The  projections  of  a  noisy  signal  are  computed  and  segmented.  The 
individual  segments  are  padded  with  zeros  and  Fourier  transformed.  The  power  spectra  of 
the  segments  of  the  projection  are  averaged  to  derive  an  estimate  of  the  power  spectrum  of 
that  one  projection.  The  same  procedure  is  carried  out  for  each  projection  to  generate  the 
2-0  power  spectrum  estimate. 

The  Blackman-Tukey  algorithm  derives  a  spectral  estimate  via  the  Wiener-Khintchine 
theorem,  i.e.  the  power  spectrum  of  a  stochastic  signal  is  the  Fourier  transform  of  its 


autocorrelation.  For  a  sampled  1-0  signal,  the  autocorrelation  is  computed  for  a  number  or 
allowed  lags  (shifts)  and  Fourier  transformed.  For  a  2-0  signal,  the  calculation  of  the  2-0 
autocorrelation  makes  this  approach  computationally  expensive.  However,  once  the 
projections  have  been  derived,  this  approach  can  be  performed  in  1-D  rapidly  and  cheaply. 
By  the  filter  theorem,  the  projection  of  the  autocorrelation  is  the  autocorrelation  of  the 
projections.  The  1-D  autocorrelation  of  each  projection  can  be  rapidly  computed,  Fourier 
transformed,  and  displayed  in  sinogram  or  polar  format  to  give  an  estimate  of  the  2-D 
power  spectrum. 

II.A.6.  Linear,  Space-Variant  Operations 

In  recent  years,  a  considerable  amount  of  effort  has  been  directed  at  developing 
optical  methods  of  implementing  space-variant  operations,  in  order  to  broaden  the 
applicability  of  optical  processing.  For  a  review  of  this  work,  see  Goodman  ^9B1).  It  is 
natural,  therefore,  for  us  to  investigate  the  application  of  the  Radon  transform  to  such 
operations.  <Ve  will  see  that  Radon-space  implementation  of  general  space-variant 
operations,  though  theoretically  possible,  usually  offers  little  if  any  advantage  over  direct 
processing.  For  some  special  cases,  however,  the  Radon  approach  can  be  very  useful. 

A  general  linear,  space-variant  operation  on  a  2-D  function  f(r)  may  be  expressed  as  a 
superposition  integral: 


r  r 

;t)  =  d*r'  f(r')  h(r;r'). 


where  the  kernel  h(r;r')  can  be  regarded  as  a  space-variant  impulse  response.  Since  the 
superposition  kernel  is  a  function  of  both  the  input  and  output  coordinates,  and  is  therefore 
4-D,  we  cannot  derive  unique  1-D  projections  of  h( r;r’ )  in  the  manner  described  by  eq.  ,2;. 
.Ve  could  derive  a  generalized  proiection  X^p.^jp',*')  of  h(r;r')  by  integration  over  the 
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.nput  and  output  variables  and  examine  the  relationship  between  ;r  and  that  fields  '< 

>Ve  have  already  seen  some  cases,  e.g.  the  Fourier  and  Hartley  transforms,  wnere  the  close 
kinship  or  the  space-variant  integral  kernel  and  the  Radon  kernel  allow  the  operations  to  oe 
directly  performed  in  this  manner.  But  for  the  general  space-variant  operation,  we  will 
instead  consider  an  alternative  treatment  made  by  Bamier  and  Hofer-Alfeis  v14d2).  They 
proved  that  2-D  space-variant  operations  can  be  considered  to  be  a  special  case  of  4-0 
space-invariant  convolution,  i.e. 

3(r)  =  g'(r;r'=0)  =  [f'(r;r')  h(r;r')J  |r._0  ,  ,31 ) 

where  f'(r;r')  =  f(r)  j(r  +  r'),  and  the  operator  denotes  4-D  convolution.  Deriving 
f'(r;r')  involves  sampling  a  4-0  smeared  version  of  f(r),  and  so  is  somewhat  akin  to  back- 
proiection.  Bamier  and  Hofer-Alfeis  proposed  a  means  of  implementing  the  4-D 
convolution  optically  via  sequential  2-D  convolutions  for  the  case  of  a  bandlimited  space- 
variant  impulse  response.  By  extension  of  the  filter  theorem  [eq.  t13)]  to  4-D,  the 
convolution  can  theoretically  be  performed  via  1-D  convolutions  in  Radon  space  once  the 
projections  of  the  4-D  functions  have  been  derived.  The  4-D  generalization  of  the 
projection  operation  [eq.  (2)]  is  obtained  in  analogous  fashion  to  the  3-D  case  Section 
I.B.I.),  i.e.  the  1-D  projection  of  a  4-D  function  is  generated  by  integration  over  the  3-D 
volume  normal  to  the  4-D  unit  vector  defining  the  azimuth  of  the  protection.  Three  angles 
,3,3,^/  are  required  to  specify  this  unit  normal.  For  clarity,  we  respecify  the  arguments 
,r;r')  of  the  4-D  functions  by  the  notation  (rj,  where  the  subscript  denotes  tne 
dimensionality  of  the  space.  Similarly,  we  define  the  4-D  volume  element  d*r  =  d1r  dlr'. 
The  1-D  protection  of  the  4-D  input  function  f‘(rt)  is  therefore: 


kf",  p,a,3,YJ 


32; 


d*r  f'{r»i  5[p  *  r,  •  n*  J  . 


Note  that  the  definition  of  f'(r„j  =  f ( r )  d(r  ♦  r‘)  allows  some  simplification  or  this 
expression  by  evaluating  the  integral  over  d*r‘.  However,  the  protection  of  the  kernel  h(r»; 
cannot  be  so  simplified,  in  general.  Extending  the  filter  theorem  [eq.  (13)]  to  4-0,  we 
have 

g'(rj  =  R*'l[Xg.(p,a,3,Y)]  =  ( Xf,(p,a,S,Y)  *  *h(P.a,B.Y)]  ,  (33) 

where  R*"1  is  the  4-0  inverse  Radon  transform.  The  desired  output  g(r)  of  the  2-0  space- 
variant  operation  is  obtained  by  evaluation  of  g'(r»)  =  g'(r;r‘)  at  the  2-D  plane  defined  by 
r'  =  0.  Since  each  1-0  convolution  influences  every  point  in  the  4-0  convolution  (and 
hence  every  point  in  the  2-0  output  plane)  via  back-projection,  there  are  no  computational 
shortcuts — only  nonessential  4-0  output.  In  Radon  terms,  mapping  the  2-D  input  function 
to  4-0  space  and  performing  a  4-0  space-invariant  convolution  avoids  the  necessity  of 
operating  on  one  proiection  of  the  2-0  input  f(r)  with  multiple  generalized  projections  of 
the  4-0  kernel  h(r;r')  to  obtain  one  projection  of  the  2-0  output  g(r).  However, 
performing  the  forward  and  inverse  Radon  transforms  of  4-0  functions  are  very  intensive 
computational  processes  which  would  require  special-purpose  hardware  if  they  are  to  be 
performed  rapidly  and  economically.  To  illustrate  the  scope  of  the  problem,  consider  that 
the  forward  transform  requires  the  calculation  of  a  volume  integral  for  each  point  in  each 
projection.  For  a  300  x  500  input  f(r),  the  general  space-variant  kernel  h ( r ; r ' )  has  500*  = 
6.25  x  lO1"  data  points.  Calculation  of  each  of  500*  projections  requires  500  volume 
integrals  over  500J  points.  The  difficulties  of  performing  the  4-0  back-projection  are 
similarly  prodigious.  4s  will  be  discussed,  Barrett  (1981)  proposed  a  hybrid  3-0  Radon- 
space  signal  processor  that  could  be  adapted  to  these  4-0  applications,  but  the  addition  of 
one  more  dimension  significantly  complicates  the  data  storage  and  manipulation 
requirements.  Hence,  performing  the  general  space-variant  operation  in  Radon  soace  via 
the  4-0  convolution  algorithm  has  no  obvious  advantage  over  direct  digital  processing  at 


this  time. 


II.A.7.  Bilinear  and  Nonlinear  Operations 

For  1-0  signals,  a  number  of  processing  algorithms  have  been  developed  that  operate 
on  the  signal  in  a  multilinear  or  nonlinear  manner  for  such  purposes  as  voice  pattern 
recognition  and  echo  deconvolution.  Examples  include  coordinate- frequency 
representations  (e.g.  sliding-window  spectrum,  Woodward  ambiguity  function,  Wigner 
distribution  function  (WDF)),  triple  correlation  (Lohmann  and  Wirnitzer,  1984),  and  the 
cepstrum  (Childers  et  al.,  1977).  The  success  of  these  algorithms  for  certain  1-0  signal¬ 
processing  tasks  has  stimulated  research  into  2-D  analogs,  but  these  are  usually 
computationally  intensive  and  hence  not  often  implemented  digitally.  In  some  cases, 
optical  processing  has  been  profitably  applied,  notably  for  coherent  optical  computation  of 
the  Wigner  distribution  function  of  2-0  data  (Bamler  and  Glunder,  19d3).  Those  operations 
based  on  Fourier  transforms  (e.g.  WDF)  or  on  nonlinear  point  processing  (e.g.  cepstrum) 
may  be  implemented  in  Radon  space.  Using  a  fast  optical  Radon  transformer  and  1  -0 
analog  or  fast  digital  processing,  the  2-0  operation  may  be  performed  profitably.  An 
example  of  such  an  operation  is  coordinate-frequency  representation  of  2-D  functions. 

A  simultaneous  representation  of  the  coordinate  and  frequency  distribution  of  the 
energy  in  a  nonstationary  signal  has  proven  useful  in  a  dumber  of  applications,  e.g.  radar 
signal  processing  (Woodward,  1953)  and  speech  processing  (Oppenheim,  1970).  Such  a 
representation  is  intended  to  give  a  picture  of  the  'local*  frequency  spectrum  of  the  signal, 
i.e.  the  frequency  content  of  the  signal  arising  from  a  particular  region  of  coordinate 
space.  Obviously,  such  a  picture  requires  twice  as  many  dimensions  in  the  representation 
space  as  in  the  signal  space.  Several  such  representations  have  been  proposed.  The  most 
direct  path  to  a  locai  soectrum  is  the  complex  spectrogram  (CS),  or  sliding-window 
spectrum,  where  a  constant  window  function  is  shifted  over  the  signal  to  specify  the  region 
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to  be  Fourier  analyzed,  i.e.  for  a  1-0  signal  f(t),  the  complex  spectrogram  Sfg  is  dermea 
as: 


Sfg  (t;v )  =  |  dt'  f(t')  g*  (t'  -  t)e'2irivt'.  i34) 

t  -a» 

This  representation  is  easily  computed  by  coherent  optics  but  the  output  is  affected  as 
much  by  the  window  g(t)  as  by  the  input  f(t).  This  potential  problem  can  be  alleviated  py 
using  a  self-windowed  representation,  such  as  the  Wigner  distribution  function  (WDF), 
which  is  commonly  defined  as: 


f  • 

Wf(t;v)  =»  |  dt‘  f(t  +  j)  f*  (t  -  y)  e‘2iri'otl 


3 


f(t  ♦  y)  f  (t  2 


-  £1 


(35) 


where  -/  ,is  the  1-0  Fourier  operator  transforming  coordinate  t‘  to  frequency  v.  This 
t1*  v 

representation  was  introduced  by  vVigner  (1932)  and  introduced  into  optics  by  Pastiaans 
(1978).  Another  closely  related  function  is  the  Woodward  ambiguity  function  (AF),  which  is 
defined  as: 


Af  - 


dt  f(t  *  j) 


*  (t 


-  t)  e-2"ivt 


if(t 

L 


;3bj 


It  is  related  to  the  WOF  through  a  double  Fourier  transform.  Several  optical  methods  for 
computing  these  representations  for  1-0  functions  have  been  introduced  (Bartelt  et  al., 
’980)  Brenner  and  lohmann,  1982)  (Eichmann  and  Dong,  1982)  (Athale  et  al.,  1982). 

Generation  or  such  representations  for  2-0  functions  presents  another  problem,  since 
the  resultant  is  a  function  of  four  variables.  Generally,  2-D  slices  of  the  4-0 


representation  are  produced,  Real  input  functions  are  assumed,  eliminating  the  need  for 
corrugating  the  shifted  function.  In  addition,  the  computation  of  the  bilinear  product 
function  is  expensive  if  done  digitally,  increasing  the  motivation  for  optical  processing.  Of 
the  representations  listed,  the  WDF  is  most  readily  computed  optically,  since  the  Fourier 
transform  of  the  product  function  is  over  the  shifted  coordinate  variable  .  Optically,  the 
product  function  is  generated  by  passing  coherent  light  twice  through  a  transparency  of  the 
signal,  either  by  reflecting  an  image  of  the  transparency  onto  itself,  by  overlaying  copies 
(Bamler  and  Clunder,  1983),  or  by  imaging  onto  a  copy  (Conner  and  Li,  1985).  The  bilinear 
product  function  is  then  Fourier  transformed  to  generate  one  slice,  Wf  (r,;p).  Shifting  the 
position  of  the  input  functions  generates  2-0  slices  for  different  values  of  r,. 

Computation  of  the  WDF  can  also  be  performed  in  Radon  space  by  talcing  projections  of 
the  optically  derived  bilinear  product  and  Fourier  transforming  in  1-0.  Easton  et  al.  (1984) 
demonstrated  generation  of  1-0  central  slices  of  the  squared  modulus  of  the  4-0  WDF  and 
later  used  the  1-0  SAW  complex  Fourier  transformer  to  produce  bipolar  2-0  slices  of  the 
4-D  WDF  of  a  2-0  real  function.  An  example  is  shown  in  Figure  21. 

II.B.  Operations  on  3-0  Signals 

Earlier,  we  stated  that  we  would  emphasize  processing  of  2-0  signals  via  a 
tomographic  transform.  However,  it  may  be  even  more  profitable  to  use  the  Radon 
transform  to  reduce  3-0  problems  to  1-0  operations,  since  digital  data  manipulation  is  even 
more  time-consuming  in  that  case.  Two  kinds  of  3-D  problems  will  be  discussed:  3-0 
purely  spatial  data,  and  2-0  spatial  data  with  a  third  dimension  (e.g.  time  or  spectrum). 
We  shall  briefly  describe  the  required  operations,  and  suggest  potential  applications. 

II.B.1.  3-0  Spatial  and  2-0  Spatial  >1-0  Temporal  Signals 

In  Section  1.8.1.  we  described  the  decomposition  of  a  3-D  function  into  a  set  or  1-0 
proiections  by  integration  over  parallel  planes.  The  protection  operation  is  identical  to  eq. 


■  2'.,  except  that  the  delta  function  reduces  the  volume  integral  to  a  set  or  planar  integrals. 
Given  a  3-D  function  f(x,y,z),  we  wish  to  integrate  the  function  over  a  set  of  planes 
normal  to  the  3-D  unit  vector  n.  Two  angles  are  required  to  define  the  normal  to  a  plane, 
commonly  the  azimuth  $  and  the  colatitude  3.  The  displacement  of  the  parallel  plane  from 
the  origin  is  again  defined  as  p,  so  that  n  =  p/p.  The  3-D  projection  operation  can  then  be 
expressed: 


x  p =  1 


r  «  r  ®  t  m 


I  I  I 

J  -oo>  —  <m 


d’r  f(r)  5(p  -  r*n)  . 
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The  3-D  version  of  the  central-slice  theorem  states  that  the  1-D  Fourier  transform  of  a 
proiection  of  a  3-D  function  yields  one  line  through  the  origin  of  the  3-D  Fourier 
transform. 

The  3-D  back-projection  operation  is  again  very  similar  to  the  2-D  case,  but  now  the 
1-D  function  is  smeared  over  the  original  proiection  plane  normal  to  n.  Kepeating  this  for 
all  directions  n  generates  a  3-D  summation  image  b(r).  In  the  filtering  step,  however, 
there  is  a  significant  qualitative  difference  between  the  2-D  and  the  3-D  cases.  Recall 
that  in  2-D,  the  Fourier  space  filter  for  the  1-D  projection  is  H(v)  =  |v|,  and  the 

Coordinate  space  counterpart  is  h(p)  =  -r— j— j,  which  falls  off  slowly  with  p.  The 

2lT  p 

corresponding  filter  for  the  inverse  3-D  Radon  transform  is  H(<j)  =  2ia2,  where  j  is  the 
magnitude  of  the  3-D  spatial  frequency  vector  (£,n,0  (Barrett,  1981).  The  coordinate 
space  filter  is  easily  found,  since  multiplication  by  -(2ir2p2)  in  the  frequency  domain 
corresponds  to  taking  the  laplacian  in  the  space  domain  (Gaskill,  1978).  The  expression  for 
the  inverse  3-D  Radon  transform  is  therefore: 


f(r)  =  *27  7iWr>] 


38) 


where  b(r)  is  the  3-D  summation  image.  Filtering  for  the  3-D  inverse  transform  is 


tnerefore  a  local  operation,  in  contrast  with  the  2-0  case. 

Barrett  vl98V,  proposed  a  hybrid  3-0  Radon-space  signal  processor  composed  or  an 
optical  system  to  derive  the  projections,  digital  storage,  1-0  signal  processing,  and  optical 
back-proiection.  The  input  function  was  assumed  to  be  a  collection  of  2-D  image  frames, 
i.e.  a  movie  film,  where  each  frame  is  assumed  to  be  a  ‘slice'  through  the  3-0  obiect.  All 
of  the  3-0  versions  of  the  operations  described  in  section  II. A.  could  be  performed  by  this 
system,  including  3-0  Fourier  transformation  and  convolution.  Such  a  system  should  be 
capaole  of  performing  3-0  complex  Fourier  transforms  on  5G0J  data  points  in  less  than  4 
hours.  A  digital  system  common  at  the  time  (POP  11/34  +  array  processor)  would  have 
required  two  days. 

Such  a  system  can  also  be  applied  to  2-0  spatial  ■►1-0  temporal  signals  (e.g.  movies) 
for  joint  spatial/temporal  filtering.  A  possible  application  would  be  to  stellar  speckle 
mterrerometry,  allowing  the  averaging  filter  impulse  response  to  vary  temporally.  Such 
operations  are  feasible  by  digital  means,  but  are  expensive  and  time-consuming. 

II.B.2.  2-0  Spatial  ♦  1-0  Spectral  Data 

Optical  detection  and  display  systems  are  best-suited  to  2-0  data  formats.  In  white- 
light  images,  a  third  dimension  of  information  has  been  encoded  in  the  spectrum  of  each 
image  point.  The  Radon  transform  provides  a  mechanism  by  which  we  may  use  2-0 
detectors,  signal  processors,  and  display  devices  to  manipulate  the  spectral  data  while 
retaining  the  ability  to  regenerate  the  image.  For  example,  if  we  have  a  white  light  2-0 
image,  we  can  derive  the  set  of  1-0  projections  of  that  image  as  described  in  section 
I.C.1.  The  1-0  projections  can  be  spectrally  dispersed  in  the  orthogonal  dimension, 
allowing  2-0  filtering  to  be  performed  on  the  joint  spatial/spectral  protection.  The  2-0 
filtered  signal  can  be  'inversely  dispersed",  to  rederive  spectrally-filtered  1-0  protections, 
and  a  2-D  filtered  image  then  reconstructed  via  the  inverse  Radon  transform,  Sucn  a 


svstem  could  oe  used  ror  spectrai  matched  riltered  imaging  JLohmann  and  Maul,  1981)  ,Yu, 
1984)  or  imaging  spectroscopy. 

III.  Summary  and  Conclusions 

We  have  discussed  the  reduction  of  2-D  signal  processing  operations  to  1-0  operations 
via  the  Radon  transform  for  the  purpose  of  gaining  flexibility,  precision,  and  mechanical 
advantages  over  direct  optical  signal  processing.  This  technique  is  most  readily  applicable 
to  operations  based  on  Fourier  transforms  and  convolution.  Several  optical  systems  were 
discussed  that  are  capable  of  performing  the  forward  and  inverse  dimensional 
transformations,  and  a  number  of  applications  were  considered,  some  already  demonstrated 
and  some  postulated.  The  authors  believe  that  many  of  the  fruits  of  this  technique  have 
yet  to  be  harvested,  and  we  encourage  workers  in  signal  processing  to  investigate  the 
utility  of  protection  operations  in  their  own  applications. 
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Two-dimensional  complex  Fourier  transform  via  the 
Radon  transform 

Roger  L.  Easton,  Jr.,  A.  J.  Ticknor,  and  H.  H.  Barrett 


A  hybrid  system  has  been  constructed  to  perform  the  complex  Fourier  transform  of  real  2  D  data.  The 
system  is  based  on  the  Radon  transform;  i.e..  operations  are  performed  on  1-D  projections  of  the  data  The 
projections  are  derived  optically  from  transmissive  or  reflective  objects,  and  the  complex  Fourier  transform  is 
performed  with  SAW  filters  via  the  chirp  transform  algorithm.  The  real  and  imaginary  parts  of  the  2-D 
transform  are  produced  in  two  bipolar  output  channels. 


through  the  2-D  transform.  When  the  transforms  are 
plotted  in  polar  format  on  CRT  screens,  the  real  and 
imagination  parts  of  the  2-D  Fourier  transform  of  the 
object  are  displayed.  The  system  can  be  used  with 
either  transmissive  or  reflective  input  data,  and  the 
illumination  may  be  incoherent.  We  have  previously 
reported  on  the  application  of  this  system  to  computa¬ 
tion  of  power  spectra,10  but  results  of  2-D  complex 
transformation  are  given  here  for  the  first  time. 

11.  Optical  Fourier  Transformation 

Fourier  transformation  by  coherent  optics  has  been 
the  basis  of  optical  processing  for  many  years  and 
found  use  even  before  invention  of  the  laser.  Coheren  t 
optica]  systems  can  compute  the  squared  modulus  of 
the  Fourier  transform  virtually  instantaneously  but 
are  limited  in  performance  by  speckle  noise  and  by  the 
available  input  transducers  (spatial  light  modulators). 
Using  the  proper  optical  configuration,  it  is  easy  to 
show  that  the  correct  amplitude  and  phase  of  the 
transform  are  produced  at  the  output  plane  <  limited  by- 
aberrations  in  the  transform  lens),  but  the  necessity  of 
square-law  detection  makes  separation  of  the  ampli¬ 
tude  and  phase  components  of  the  transform  lor.  near¬ 
ly  equivalently,  of  the  real  and  imaginary  parts)  diffi¬ 
cult. 

A  considerable  body  of  work  has  been  done  on  pro¬ 
duction  of  the  Fourier  transform  by  incoherent  optics 
with  the  aim  of  gaining  significant  advantages  over 
coherent  optics  in  output  noise  and  flexibility  of  inputs 
while  retaining  the  speed  advantage  over  digital  com¬ 
putation.  Katvl1  used  a  temporally  incoherent  source 
in  the  coherent  optics  format  with  appropriate  disper¬ 
sion  correction.  The  requirement  for  spatial  coher¬ 
ence  remains,  and  derivation  of  the  complex  transform 
is  difficult.  Other  systems  use  geometric  shadow¬ 
casting  to  image  the  input  on  a  reference  mask  of 
known  spatial  frequency  and  phase.  The  integrated 
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I.  Introduction 

The  utility  of  the  2-D  Fourier  transform  as  a  tool  for 
signal  processing  is  well  known.  Its  computation  is 
usually  performed  digitally  or  by  coherent  optics. 
Other  techniques  have  been  demonstrated  to  compute 
the  2-D  transform  using  incoherent  illumination.1'9 
Each  of  these  methods  has  inherent  advantages  and 
disadvantages.  Digital  computation  on  a  general  pur¬ 
pose  computer  is  precise  but  slow,  even  with  the  FFT 
algorithm.  In  addition,  it  can  suffer  from  aliasing 
problems  if  the  data  are  inadequately  sampled.  The 
use  of  special  purpose  hardware,  3uch  as  array  proces¬ 
sors,  can  speed  the  process  considerably,  but  digital 
techniques  cannot  as  yet  approach  video  rates  (30 
frames/sec)  with  large  arrays.  Optical  methods  to 
compute  the  Fourier  transform  have  been  developed, 
but  each  has  disadvantages  limiting  its  utility.  We 
have  constructed  a  system  capable  of  performing  com¬ 
plex  Fourier  transforms  of  2-D  input  data  at  video 
rates.  The  system  is  based  on  the  Radon  transform 
and  the  chirp  Fourier  transform.  An  optical  scanner 
produces  1-D  projections  of  the  input  data,  which  are 
Fourier  transformed  in  1-D  by  a  surface  acoustic  wave 
chirp  transformer.  The  real  and  imaginary  parts  of 
the  transform  are  produced  simultaneously  in  separate 
channels.  A  single  projection  is  derived  in  10  ^sec,  and 
the  complex  transform  is  produced  <30  usee  after 
commencement  of  the  scan.  By  the  central  slice  theo¬ 
rem.  these  1-D  transforms  are  equivalent  to  lines 
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light  at  one  position  in  the  output  plane  is  proportional 
to  the  Fourier  coefficient  at  one  spatial  frequency. 
The  technique  of  Mertz,3  later  refined  by  Richardson,  5 
produced  the  reference  masks  via  the  moire  pattern 
created  by  two  Fresnel  zone  plates.  By  sequential 
replacement  of  the  second  zone  plate  with  one  with 
spatial  frequencies  in  quadrature.  Richardson  was  able 
to  compute  the  cosine  and  sine  transforms  separately 
giving  the  real  and  imaginary  parts  of  the  Fourier 
transform.  This  implementation  can  be  analyzed  as 
the  chirp  algorithm  for  Fourier  transformation,  which 
decomposes  the  correlation  with  the  Fourier  kernel 
into  multiplication  and  convolution  with  quadratic 
phase  factors  as  described  later  in  this  paper.  Leifer  et 
al.*  used  a  stored  reference  mask  with  a  limited  range 
of  spatial  frequency  and  orientation  in  a  shadow-cast¬ 
ing  correlator  for  alphabetic  character  recognition.  In 
all  these  systems,  the  shift  required  for  correlation  of 
the  input  with  the  reference  spatial  frequency  mask  is 
accomplished  by  optical  parallax,  and  no  physical 
movement  is  required.  The  maximum  spatial  fre¬ 
quency  response  of  these  systems  is  limited  by  vignett¬ 
ing  of  the  reference  masks  and  by  diffraction  (since 
geometric  optics  is  assumed).  The  vignetting  problem 
may  be  solved  by  using  a  moving  correlator  at  the 
expense  of  slower  calculation  and  increased  complex¬ 
ity.  Even  here,  the  scanning  need  not  be  physical 
motion  if  an  imaging  detector  is  used.5  However,  the 
geometric  optics  assumption  severely  limits  the  spatia* 
frequency  response  of  these  incoherent  correlation 
systems  to  arrays  of  100  X  100  pixels  or  so.  In  addi¬ 
tion,  the  spurious  terms  present  in  the  output  plane 
decrease  contrast  and  reduce  output  dynamic  range. 

Other  authors  have  investigated  different  avenues 
to  Fourier  transform  computation.  Recent  work  by 
Tai  and  Aleksoff’5  has  demonstrated  production  of 
complex  transforms  of  incoherently  illuminated  data 
by  selection  of  the  proper  output  term  from  a  grating 
interferometer.  This  approach  is  limited  to  1-D  data, 
however.  Xu  et  al~  have  produced  the  complex  trans¬ 
form  of  incoherently  illuminated  2-D  data  occupying 
one-half  of  the  input  plane.  A  symmetric  object  is 
synthesized  by  reflection  through  the  origin  and  pro¬ 
cessed  through  two  illumination  channels  polarized 
orthogonally.  The  system  performs  well,  but  the  re¬ 
striction  on  input  format  limits  its  utility.  George  and 
Wang’ also  have  performed  Fourier  cosine  transforma¬ 
tion  of  transmissive  or  reflective  objects  in  incoherent 
light  by  synthesis  of  a  symmetric  object  followed  by  an 
achromatic  optical  Fourier  transform.  A  double  im¬ 
age  of  the  input  is  produced  interferometrically,  and 
the  output  of  the  optical  system  is  the  cosine  transform 
on  a  bias.  Adjustment  of  the  interferometer  allows 
separate  generation  of  the  sine  transform.  The  output 
signal  is  detected  with  a  photodiode  array  for  later 
digital  manipulation.  The  bias  could  be  subtracted 
electronically  or  interferometrically.  They  report  sys¬ 
tem  response  to  20  cycles/ mm.  and  their  results  agree 
very  well  with  calculations.  This  svstem  has  the  po¬ 
tential  disadvantage  of  nonsimultaneous  generation  of 
the  cosine  and  sine  transforms. 
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Glaser  et  al.9  have  implemented  the  chirp  transform 
algorithm  optically  to  produce  the  complex  transform. 
A  holographic  filter  is  used  to  perform  the  convolution 
with  the  quadratic  phase  factor,  thus  requiring  a  tem¬ 
porally  quasi-coherent  source.  Spatial  coherence  is 
not  required.  The  optical  output  is  a  spatial  carrier 
modulated  by  the  complex  transform,  from  which  the 
real  and  imaginary  parts  may,  in  principle,  be  derived 
by  digital  demodulation  at  the  cost  of  temporal  pro¬ 
cessing  capacity. 

III.  Radon-Fourier  Transformer 

All  the  2-D  Fourier  transforming  systems  discussed 
above  are  restricted  in  utility  by  limitations  on  speed, 
format  of  input  and/or  output,  space-bandwidth  prod  - 
uct,  or  dynamic  range.  Many  of  these  systems  have 
proven  useful  in  some  applications,  but  none  truly  fills 
the  need  for  rapid  calculation  of  the  complex  2-D  Fou¬ 
rier  transform  with  large  space-bandwidth  product. 
Using  a  different  principle,  we  have  constructed  a  sys¬ 
tem  which  can  potentially  compute  complex  Fourier 
transforms  of  large  arrays  at  video  rates.  The  complex 
transform  is  generated  as  cosine  and  sine  transforms, 
i.e.,  the  real  and  imaginary  parts  of  the  transform. 
The  two  outputs  are  obtained  simultaneously.  Opera¬ 
tion  is  based  on  the  Radon  transform,11  which  decom¬ 
poses  a  function  of  Af-dimensions  into  the  complete  set 
of  1-D  projections  by  integration  over  .V/  —  1  dimen¬ 
sions.  For  the  2-D  case,  projections  are  obtained  by- 
integration  over  sets  of  parallel  lines.  The  primary- 
theorem  of  the  Radon  transform  states  that  a  function 
can  be  reconstructed  from  the  complete  set  of  its  pro¬ 
jections  and  serves  as  the  operating  principle  of  medi¬ 
cal  computed  tomography.  The  Radon  transform  has 
also  been  shown  to  be  useful  in  general  signal  process¬ 
ing,  including  pattern  recognition, 13  image  filter¬ 
ing, U1S  bandwidth  compression.16-1"  computation  of 
the  Wigner  distribution  function.16  and  Fourier  spec¬ 
trum  analysis.10  19--’° 

The  utility  of  the  Radon  transform  for  signal  pro¬ 
cessing  is  due  to  the  central-slice,  or  projection-slice, 
theorem,  which  states  that  the  1-D  Fourier  transform 
of  a  1-D  Radon  projection  yields  one  line  through  the 
2-D  Fourier  transform  of  the  2-D  function.  The  1-D 
transform  passes  through  the  origin  of  2-D  Fourier 
space,  and  its  orientation  is  determined  by  the  orienta¬ 
tion  of  the  lines  of  integration.  Since  systems  exist 
that  can  rapidly  compute  1-D  Fourier  transforms  i  e.g.. 
CCD,  SAW.  or  AO),  adopting  the  Radon  transform 
approach  makes  possible  rapid  computation  of  the  2-D 
Fourier  transform. 

The  system  for  producing  the  Radon  transform  of 
the  2-D  data  has  been  discussed  previously.- 
Suffice  it  to  say  that  the  projections  of  fhe  2-D  distri¬ 
bution  of  intensity-  transmission  t of  a  transparency  i  or 
reflectance  ifor  reflective  objects i  ..re  derived  by  pro¬ 
jecting  a  line  of  light  on  the  input  plane  and  integrating 
the  light  transmitted  or  reflected  with  a  detector  The 
output  of  the  detector  is  proportional  to  the  line  inte¬ 
gral  of  transmission  or  reflectance.  Sweeping  the  line 
of  light  perpendicular  to  itself  across  the  input  data 
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produces  a  temporal  signal  proportional  to  one  line- 
integral  projection.  Rotation  of  the  azimuth  of  sweep 
with  a  prism  allows  production  of  the  complete  Radon 
transform  as  a  sequence  of  1-D  temporal  signals  out  of 
the  detector.  For  obvious  reasons,  this  optical  system 
is  termed  a  flying-line  scanner. 

From  the  central-slice  theorem  mentioned  above, 
the  1-D  transform  of  a  projection  is  one  line  through  a 
polar  plot  of  the  2-D  transform.  Farhat  et  al.2X  have 
adapted  both  coherent  optical  transformation  and 
shadow-casting  correlation  to  perform  the  1-D  Fourier 
transformation  of  Radon- transformed  data.  Their  in¬ 
coherent  transformer  produces  full  complex  trans¬ 
forms  of  complex  input  data  by  using  two-color  chan¬ 
nels.  We  have  taken  a  different  tack,  disposing  of 
optical  Fourier  transformation  altogether,  and  instead 
implementing  the  chirp  Fourier  transform  algorithm 
with  surface  acoustic  wave  filters.  The  chirp  trans¬ 
form  results  from  a  decomposition  of  the  Fourier  ker¬ 
nel: 

exp(-2rii'£)  »  jexp(-ir^j^  j  X  lexp[— ix(/3t)2! 


:  ^P  [‘*(j -*)]} 


Thus  the  Fourier  transform  may  be  written 


HD1 


fit)  exp<— 2 rmt)dt  •  exp 


X  j  {/(f)  exp[— ir(/Jt)2][  X  expj^r  -  dtj  dt 
»  (expHrUit)*] 

X  (i/(t)  exp(— irf/Jt)2]!  .  exp[ir(/3t)ll)I|(^tfa,),  (2) 

where  *  denotes  convolution.  The  complex  exponen¬ 
tials  are  linear  FM  signals,  i.e.,  the  frequency  varies 
linearly  with  time,  and  have  been  named  chirps  by  the 
radar  community.  They  are  also  called  quadratic 
phase  factors  for  obvious  reasons.  The  instantaneous 
frequency  of  the  positive  complex  quadratic  phase 
term  exp[+i7r(/Jt)2]  at  time  t„  is 

-Lx**lu  -  + 

2jt  dt 2 

which  increases  with  tn.  Hence,  it  is  called  an  upchirp, 
while  the  negative  exponential  is  a  downshirp.  Using 
the  three  chirps,  Fourier  transformation  can  be  broken 
down  into  the  following  3teps: 

(1)  multiplication  of  the  signal  by  exp(— iird2t2),  a 
downchirp; 

(2)  convolution  of  the  product  with  an  upchirp, 
expi+iir&t2); 

(3)  multiplication  of  the  filtered  signal  by  a  down- 
chirp,  exp(-iird2t2). 

The  parameter  r d2  in  the  chirp  signal  is  called  the 
chirp  rate  and  is  the  same  for  all  three  chirps.  The 
temporal  output  signal  is  a  scaled  version  of  the  Fouri¬ 
er  transform  of  the  temporal  input  signal  with  the 


temporal  frequency  related  to  the  output  temporal 
coordinate  via  u  =  &2t.  If  only  the  modulus  of  the 
transform  is  required,  the  third  step  can  be  deleted. 
This  analysis  assumes  that  the  chirps  are  complex  and 
of  infinite  length. 

One  way  to  obtain  the  chirp  impulse  responses  nec¬ 
essary  to  implement  the  transform  is  via  surface  acous¬ 
tic  wave  (SAW)  chirp  filters.22-24  A  SAW  filter  con¬ 
sists  of  a  piezoelectric  crystalline  substrate  on  which 
two  aluminum  interdigital  transducers  (IDTs)  have 
been  deposited.  When  a  signal  is  applied  to  the  input 
IDT,  the  electric  field  across  transducer  fingers  of  op¬ 
posite  polarity  generates  a  deformation  of  the  crystal 
surface  via  the  piezoelectric  effect.  The  deformation 
travels  along  the  crystal  surface  as  a  sound  wave.  At 
the  output  IDT,  an  electronic  signal  is  regenerated 
from  the  sound  wave  by  the  inverse  piezoelectric  inter¬ 
action.  By  proper  design  of  the  separations  and  over¬ 
laps  of  the  fingers  in  the  IDTs,  any  of  a  wide  variety  of 
impulse  responses  can  be  generated.  For  a  chirp  filter, 
the  separations  of  the  fingers  are  varied  to  obtain  an 
impulse  response  h(t)  whose  frequency  increases  or 
decreases  from  some  initial  carrier  frequency  wo  at  rate 
a,  i.e.,  of  the  form 


hit)  «■  cos 


[Ht)]' 


Again,  the  instantaneous  frequency  of  the  chirp  at 
time  t„  is 


wojfcat, 

2r 


,aCn 


As  before,  the  positive  term  is  called  an  unchirp.  Real¬ 
istically,  the  filter  must  have  a  finite  temporal  re¬ 
sponse,  so  the  cosine  function  must  be  windowed  by  a 
function  with  compact  support,  e.g.,  a  Rect  function  or 
a  Hamming  window.  For  chirp  Fourier  transforma¬ 
tion,  the  premultiplication  and  postmultiplication 
chirps  can  be  generated  by  applying  an  impulse  input 
to  SAW  chirp  filters  of  the  proper  sign  (i.e.,  upchirp  or 
downchirp).  The  convolution  is  performed  by  apply¬ 
ing  the  signal  to  a  similar  filter.  It  is  important  to  note 
that  the  SAW  chirp  filter  impulse  response  is  a  real 
function  of  the  form  Ait)  cos(w£  ±  at2),  not  the  com¬ 
plex  exponential  seen  above.  The  function  Ait)  is  an 
apodization  of  the  chirp,  necessitated  by  the  finite 
output  signal  length,  and  w  is  the  initial  angular  fre¬ 
quency  of  the  chirp.  The  chirp  transform  algorithm 
may  still  be  implemented,23  but  the  steps  now  become 

(1)  premultiplication  by  a  downchirp 

Rect  x  cos  t  -  • 

(2)  filter  this  signal  with  impulse  response 

Alt)  =*  Rect  — —  ^  X  cos^„  t  +  ' 

(3)  postmultiply  the  filtered  signal  by  two  up- 
chirps  separately: 
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(a)  Rect  ^  x  cos (^i+t  + 

lb)  Rect  ^ - — —  -  ^  x  cos  ^ 

*  Rect  X  sin  ^>+t  +  • 

(4)  low-pass  filter  both  outputs  of  step  (3).  The 
terms  r+  and  r_  represent  the  temporal  width  of  the 
upchirp  and  downchirp,  respectively,  and  are  also 
called  the  time  dispersions  of  the  chirps.  Similarly,  o>+ 
and  <j>-  are  the  initial  angular  frequency  of  the  upchirp 
and  downchirp. 

The  output  temporal  frequency  is  related  to  the 
temporal  position  in  the  output  signal  by  the  relation 

v  m  (o>+  —  +  at)/2r.  (3) 

Since  the  carrier  frequencies  in  SAW  filters  are  in 
the  rf  regime  (w  s *  15-300  MHz),  multiplications  can 
be  performed  in  rf  mixers.  The  discrepancy  in  the  sign 
of  the  two  postmultiplication  SAW  chirps  relative  to 
that  in  the  complex  chirp  algorithm  given  above  results 
from  rf  double-sideband  mixer  multiplication.  Such 
mixers  yield  product  terms  as  modulations  on  carriers 
at  the  sum  and  difference  frequencies  of  the  original 
carriers.  That  is,  given  two  signals  A(t)  and  B(t)  mod¬ 
ulating  carriers  at  angular  frequencies  <*i0  and  <*>(,,  re¬ 
spectively,  the  action  of  the  rf  mixer  is  to  produce  an 
output: 

-4(f)  C09<uia  f)  X  B(t )  coslwj  £)  • 

X  |coa[(w0  +  Wjjt]  +  cos((uia  -  UI6)C]|.  (4) 

The  low- pass  filter  selects  the  difference  frequency 
term,  and  so  the  sign  of  the  postmultiplication  chirp 
must  be  the  same  as  that  of  the  convolution  filter  to 
obtain  demodulation.  The  signal  postmultiplied  by 
the  cosine  upchirp  is  the  real  part  of  the  transform, 
while  that  multiplied  by  the  sine  upchirp  is  the  imagi¬ 
nary  part  of  the  transform. 

For  maximum  time-bandwidth  product  in  the  out¬ 
put  signal,  the  requirements  on  the  chirps  are  that  the 
time  dispersions  of  the  premultiplication  and  the  con¬ 
volution  chirp  be  related  by  r_  ■  r+/2,  and  that  the 
bandwidth  of  the  convolution  chirp  be  twice  that  of  the 
multiplication  chirps.23  The  two  outputs  are  propor¬ 
tional,  respectively,  to  the  real  and  imaginary  parts  of 
the  Fourier  transform  within  a  time  window  (r+/2  <  t 
<  r+).  The  corresponding  spectral  window  spans 
temporal  frequencies 


The  rectangular  Finite-length  window  of  the  convolu¬ 
tion  filter  has  the  effect  of  convolving  the  spectral 
components  with  a  sine  function,  which  limits  the 
number  of  resolvable  frequencies  in  the  spectrum  to 
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Fig.  1.  Schematic  of  the  1-D  SAW  complex  Fourier  transformer. 
The  temporal  signal  /(£)  from  the  photomultiplier  in  the  flying- line 
scanner  is  proportional  to  one  projection.  The  impulse  response  of 
the  SAW  filters  is  ht(£).  The  microcomputer  controller  sends  a 
trigger  signal  to  the  digital  delay  generator,  which  in  turn  produces  a 
L-nsec  pulse  that  is  applied  to  the  downchirp  SAW  filter.  The 
resulting  impulse- response  signal  h.U)  is  multiplied  by  the  incom¬ 
ing  projection  signal  in  a  rf  mixer.  The  product  signal  is  applied  to 
the  upchirp  SAW  filter,  and  the  output  goes  to  the  signal-input  port 
of  the  rf  phase  comparator.  After  a  delay  of  14  nsec  ( l-nsec  resolu¬ 
tion),  the  digital  delay  generator  outputs  a  second  l-nsec  pulse, 
which  is  applied  to  the  postmultiplication  SAW  filter.  An  rf  phase 
shifter  at  the  SAW  filter  output  allows  fine  adjustment  of  the  post¬ 
multiplication  timing.  This  signal  is  applied  to  the  reference  port  of 
the  phase  comparator.  After  low  pass  filtering,  the  in-phase  I 
output  of  the  phase  comparator  is  proportional  to  the  real  part  of  the 
Fourier  transform  F(u)  (i.e.,  cosine  transform)  of  the  input  signal 
f(t).  Similarly,  the  output  of  the  quadrature  port  Q  of  the  phase 
comparator  is  proportional  to  the  imginary  part  of  F{ »).  ( i.e..  sine 
transform). 


one-fourth  of  the  time-bandwidth  product  of  the  con¬ 
volution  filter.23  SAW  chirp  filters  with  other  window 
functions  (e.g.,  Hamming)  are  available  if  smaller  side- 
lobes  are  desired  in  the  output  signal.  If  only  the 
squared-modulus  of  the  Fourier  transform  is  required, 
square- law  envelope  detection  can  be  substituted  for 
steps  (3)  and  (4).  This  is  the  algorithm  we  have  used 
previously  to  perform  2-D  spectrum  analysis  in  Radon 
space.10 

The  complex  transform  algorithm  was  implemented 
as  shown  in  Fig.  1  using  SAW  chirp  filters  from  Ander¬ 
sen  Laboratories  (models  DS-120-10-20-251A  and 
-252A),  which  have  bandwidths  of  10  MHz,  maximum 
time  dispersions  of  20  Msec,  and  a  resulting  time-band¬ 
width  product  of  200.  The  chirp  rate  a  =  2i  X  10 
MHz/20  Msec  =  iX  1012  Hz2.  The  filter  windows  were 
unweighted.  A  flying-line  scan,  producing  one  Radon 
projection,  is  made  in  10  ^sec  and  is  synchronized  with 
the  signal  driving  the  premultiplication  impulse  gener¬ 
ator  so  that  the  center  of  the  scan  is  mixed  with  the 
center  of  the  premultiplication  downchirp.  This 
time-gates  the  premultiplication  signal  for  a  maximum 
system  time-bandwidth  product.  After  Filtering  in 
the  upchirp  SAW.  the  signal  is  coherently  demodulat¬ 
ed  by  the  postmultiplication  upchirp.  To  obtain  40 
dB  of  rejection  of  the  signal  from  one  channel  of  the 
transform  from  the  other  channel,  the  time  of  the 
premultiplication  chirp  impulse  must  be  synchronized 
to  the  postmultiplication  impulse  to  an  accuracy  of 


better  than  100  psec.23  The  timing  interval  between 
impulse  inputs  to  the  premultiplication  and  postmul¬ 
tiplication  chirps  is  the  value  t'  in  step  (3).  A  digital 
delay  generator  is  used  to  provide  the  impulse  to  the 
postmultiplication  chirp  filter  with  temporal  resolu¬ 
tion  of  1  nsec.  More  precise  timing  is  provided  by 
shifting  the  phase  of  the  postmultiplication  chirp  with 
a  continuously  adjustable  rf  hybrid  phase  shifter.  The 
demodulation  is  accomplished  in  an  rf  phase  compara¬ 
tor,  a  four-port  device  which  produces  in-phase  and 
quadrature  mixed  signals  from  an  input  signal  and 
reference.  It  may  be  thought  of  as  a  combination  of  a 
signal  splitter,  a  quadrature  hybrid,  and  two  double¬ 
sideband  mixers.  The  filtered  signal  is  split  and 
mixed  with  in-phase  and  quadrature  components  of 


Fig.  2.  Performance  of  the  SAW  chirp  complex  Fourier  transform¬ 
er  In  each  of  the  four  cases  shown,  the  top  trace  is  the  signal  from 
the  flying-line  scanner,  i.e..  a  single  projection  of  the  2-D  input.  The 
second  and  third  traces  are  the  cosine  transform  and  sine  transform, 
respectively,  produced  by  the  SAW  chirp  transformer.  The  traces 
on  the  right-hand  side  are  a  computer  simulation  of  the  same  signal. 
The  object  was  a  grating  of  25°i  duty  cycle  in  a  circular  aperture.  In 
the  first  case,  the  grating  was  centered  in  the  aperture  creating  a 
svmmetric  signal  whose  Fourier  transform  is  purely  real.  In  the 
other  three  cases,  the  grating  was  translated  relative  to  the  circular 
aperture  giving  an  asymmetric  signal  with  a  complex  transform. 
Each  horizontal  division  in  the  oscilloscope  traces  represents  5  nsec, 
indicating  that  the  complete  transform  is  computed  within  .10  usee. 


the  postmultiplication  chirp.  After  low-pass  filtering 
in  each  channel,  the  in-phase  signal  is  the  bipolar 
cosine  transform,  and  the  quadrature  signal  is  the  bi¬ 
polar  sine  transform  (each  within  the  frequency  win¬ 
dow  and  convolved  with  the  sine  function  due  to  the 
finite  convolution  window  as  described  above). 

Using  the  SAW  filters  described,  the  chirp  trans¬ 
former  resolves  fifty  temporal  frequencies  in  the  win¬ 
dow  (|g|  <  2.5  MHz).  When  the  output  of  the  flying¬ 
line  scanner  is  applied  to  the  SAW  chirp  Fourier 
transformer,  the  spatial  frequency  scaling  depends  on 
the  scanning  3peed.  Typically,  we  scan  a  25-mm  aper¬ 
ture  in  10  Msec,  giving  a  spatial  frequency  range  of  ±  1 
cycle/mm  with  fifty  resolvable  points.  By  scanning  a 
10-mm  aperture  in  the  same  time,  the  spatial  frequen¬ 
cy  response  is  ±2.5  cycles/mm.  This  by  no  means  is 
the  limit  of  a  SAW  chirp  filter  or  optical  scanner  tech¬ 
nology.  Using  reflective-array  SAW  chirp  filters 
(RACs),  transformers  capable  of  resolving  3600  points 
within  a  60-Msec  output  window  have  been  reported.25 
Were  we  to  use  this  chirp  transformer  and  scan  a  30- 
mm  diam  aperture  in  30  Msec,  we  would  obtain  900 
resolvable  points  in  a  spatial  frequency  range  of  ±  15 
cycles/mm. 

The  performance  of  the  complex  Fourier  transform¬ 
er  for  a  1-D  signal  is  demonstrated  in  Fig.  2,  where  the 
output  is  compared  to  a  computer  simulation.  A  grat¬ 
ing  (75%  clear,  25%  opaque)  was  placed  in  a  circular 
aperture  of  20-mm  diameter  in  the  flying-line  scanner. 
The  azimuth  of  scan  was  oriented  so  that  the  line  of 
integration  was  parallel  to  the  grating  lines.  The  grat¬ 
ing  was  mounted  on  a  translation  stage  so  that  it  could 
be  shifted  within  the  circular  aperture.  Four  cases  are 
shown  for  both  the  actual  and  computed  outputs.  In 
each  example,  the  top  trace  is  the  output  of  the  flying¬ 
line  scanner,  i.e.,  the  Radon  transform  of  the  object  for 
one  azimuth.  The  second  and  third  traces  are  the 
cosine  and  sine  transform  outputs  of  the  complex  Fou¬ 
rier  transformer,  i.e.,  one  line  through  the  2-D  real  part 
and  imaginary  part,  respectively,  of  the  Fourier  trans¬ 
form  of  the  original  object,  via  the  central-slice  theo¬ 
rem.  The  scanning  time  is  10  Msec,  and  the  two  trans¬ 
forms  have  been  output  within  30  Msec  after  the 
beginning  of  the  flying- line  scan.  In  the  first  case,  the 
grating  is  centered  in  the  aperture  resulting  in  a  sym¬ 
metric  input  to  the  Fourier  transformer.  The  Fourier 
transform  of  a  symmetric  object  is  purely  real,  and 
hence  the  sine  transform  vanishes,  as  shown.  Also 
note  that  the  cosine  transform  is  bipolar  and  symmet¬ 
ric.  In  the  succeeding  three  cases,  the  grating  is  trans¬ 
lated  in  the  aperture,  resulting  in  an  asymmetric  input 
to  the  complex  transformer  and  a  nonvanishing  bipo¬ 
lar  and  antisymmetric  sine  transform.  The  actual 
transformer  output  agrees  very  well  with  the  computer 
simulations. 

To  produce  the  complex  2-D  Fourier  transform,  the 
central  slice  theorem  says  that  it  is  merely  necessary  to 
display  the  1-D  transforms  of  the  projections  in  the 
proper  polar  format.  However,  for  discrete  uniform 
sampling  along  both  the  azimuthal  and  radial  axes,  the 
Fourier  space  will  be  densely  sampled  near  the  origin 
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Fig.  3.  Two-dimensional  complex  Fourier  transforms  of  a  circular 
aperture.  The  object  was  a  single  circular  aperture  of  1.0- mm 
diameter,  as  shown  at  top.  The  letters  denote  the  origin  of  coordi¬ 
nates  (i.e.,  the  optical  axis)  for  each  case.  The  display  was  biased  up, 
so  that  zero  amplitude  is  the  brightness  level  shown  in  the  imaginary 
part  of  (A).  The  brightest  areas  represent  the  most  positive  ampli¬ 
tude  of  the  transform,  and  the  darkest  areas  represent  the  most 
negative  amplitude.  (A)  With  the  aperture  centered  at  the  origin, 
the  transform  is  purely  real.  (B),  (C)  The  aperture  was  translated 
from  the  optic  axis  by  ~1.4  and  2.4  diameters,  respectively,  produc¬ 
ing  fringes  due  to  the  constant  phase  term. 


and  sparsely  sampled  at  the  high  spatial  frequencies. 
The  function  so  obtained  is  equivalent  to  [f’iplj/lpj, 
where  p  is  the  2-D  frequency  vector  and  F(p)  is  the  2-D 
Fourier  transform  of  the  2-D  input  function  /( r)  * 
f{x,y).  The  radial  spatial  frequency  vector  p  is  always 
non-negative,  but  we  can  also  consider  a  radial  fre¬ 
quency  vector  u,  which  is  bipolar.  To  counter  the  1/1  pj 
weighting,  it  is  necessary  to  multiply  the  Fourier  trans¬ 
former  output  by  |u|  before  display.  This  V-shaped 
function  is  produced  electronically  by  passing  a  ramp 
function  through  an  absolute -value  amplifier.  Since 
the  SNR  of  the  Fourier  transformer  generally  de¬ 
creases  with  increasing  frequency,  the  V  function  is 
rolled  off  by  current-limiting  the  output  of  the  V  gen¬ 
erator.  For  a  signal  on  an  rf  carrier,  (i.e.,  the  generated 
magnitude  of  the  Fourier  transform  that  is  output  by 
the  convolution  chirp  Filter),  the  multiplication  by  |  w| 
is  easily  done  in  a  rf  mixer.10  Multiplication  of  the 
coherently  demodulated  signal  is  somewhat  more  dif¬ 
ficult  in  the  frequency  range  of  interest  (up  to  2.5 
MHz),  which  is  higher  than  most  analog  multiplier 
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Fig.  4.  Two-dimensional  complex  Fourier  transforms  of  a  grating 
in  a  circular  aperture.  The  spatial  frequency  of  the  grating  was  1.5 
cycles/mm,  with  a  duty  cycle  of  80%  and  the  aperture  diameter  was  6 
mm.  (A)  Cosine  transform  with  the  object  centered  on  the  optical 
axis  as  shown.  The  transform  is  even,  and  the  Airy  patterns  of  the 
circular  aperture  at  the  ±  1  orders  of  the  grating  are  clearly  seen.  ( B) 
Cosine  transform  after  translating  the  object  by  one-half  of  a  grating 
cycle.  The  linear  phase  term  resulting  from  the  translation  has 
inverted  the  phase  of  the  Airy  patterns.  (C),  (D)  Sine  transform  of 
the  object  after  translation  by  ±one-fourth  of  a  grating  cycle,  respec¬ 
tively,  relative  to  (A).  The  transforms  are  odd.  and  the  Airy  pat¬ 
terns  at  the  ±1  orders  are  out  of  phase.  (E),  (F)  The  aperture 
diameter  was  reduced  to  2.5  mm,  and  the  center  was  translated 
relative  the  optic  axis  by  a  sufficient  distance  (2  mm)  so  that  several 
cycles  of  the  linear  phase  are  visible  within  the  central  disk  of  the 
Airy  pattern.  (E)  is  the  real  part  of  the  transform  and  is  even.  (F>  is 
the  imaginary  part  and  is  odd.  Note  that  the  translation  was  in 
different  directions  in  the  two  cases,  so  that  the  fringe  direction 
differs. 


modules  can  handle  and  not  high  enough  for  rf  mixers. 
In  their  stead,  we  employed  the  Motorola  balanced 
modulator-demodulator  integrated  circuit  (MC1496) 
to  multiply  the  transformer  output  by  |  v  | .  The  bipolar 
signal  can  then  be  applied  to  the  z  axis  of  a  CRT  in  one 
of  two  ways:  the  signal  can  be  thresholded  at  ground 
to  display  the  complex  transformation  in  four  parts 
(positive  and  negative  real  and  positive  and  negative 
imaginary),  or  the  bipolar  signal  can  be  biased  up  to 
display  the  complete  real  or  imaginary  transform  at 
one  time.  Since  the  1-D  cosine  and  sine  transform 
signals  are  available  simultaneously,  they  can  be  dis¬ 
played  simultaneously  on  separate  CRTs. 

To  display  the  transform  in  the  polar  format,  we 
have  used  the  same  system  reported  previously. 10  The 
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Fig.  5.  Two-dimensional  complex  Fourier  transforms  of  two  circu¬ 
lar  apertures.  The  diameter  of  the  apertures  was  1  mm  with  their 
centers  separated  by  S  mm,  as  shown  at  top.  Again  the  letters 
denote  the  position  of  the  optical  axis  in  each  case.  (A)  Cosine 
transform  with  the  optic  axis  centered  on  the  object’s  axis  of  symme¬ 
try.  Note  the  phase  change  as  a  fringe  pasaes-from  the  central  lobe 
of  the  Airy  disk  to  the  first  ring.  The  faint  fringes  in  the  imaginary 
part  of  the  transform  are  due  to  wobble  in  the  image  rotating  prism. 
(B)  The  optic  axis  was  located  l  mm  above  the  symmetry  axis 
producing  fringes  perpendicular  to  those  from  the  double  aperture. 
The  cosine  transform  is  even,  and  the  sine  transform  is  odd.  (C)  The 
optic  axis  was  located  on  the  symmetry  axis  but  displaced  from  the 
center  of  symmetry  by  l  mm  multiplying  the  fringes  by  a  linear 
phase  term  of  lower  frequency. 


scan  azimuth  is  rotated  by  an  image-rotation  prism  via 
a  stepper  motor  resolving  200  steps.  The  maximum 
angular  resolution  in  the  transform  is  t/100  rad.  A 
bipolar  ramp  function  is  generated  and  weighted  in 
two  channels  by  the  sine  and  cosine  value  of  the  azi¬ 
muth  angle  of  the  3can.  The  resulting  outputs  are 
applied  to  the  x  and  y  deflections  of  the  CRT  scanning 
spot,  which  produces  a  line  scan  across  the  screen  at 
the  appropriate  angle.  After  completion  of  the  scan, 
values  of  the  sine  and  cosine  of  the  new  angle  are  read 
out  of  a  lookup  table  for  the  next  scan.  The  scanning 
spot  is  timed  to  reach  the  center  of  the  screen  when  the 
zero-frequency  output  of  the  Fourier  transformer  is 
applied  to  the  z  axis  of  the  CRT.  The  complete  2-D 
transform  can  now  be  generated  in  ~0.1-sec,  limited  by 
the  rotation  rate  of  the  stepper  motor.  To  allow  trans¬ 
formation  at  video  rates,  the  azimuth  of  the  Radon 


Fig.  6.  T wo-dimensional  complex  transforms  of  a  reflective  object. 
A  beam  splitter  was  introduced  into  the  flying-line  scanner  to  direct 
the  reflected  line  onto  the  detector.  Fourier  transformation  and 
display  were  performed  as  before.  The  object  was  a  grating  in  a 
circular  aperture  of  6-mm  diameter,  as  in  Fig.  4.  The  main  features 
of  the  transform  are  easily  seen,  i.e.,  the  location  and  phase  of  the 
Airy  patterns  on  the  first  orders  of  the  grating  spectrum.  The 
signal-to-noise  is  less  than  in  the  transmissive  case  due  to  the  lower 
reflectance  and  lower  modulation  in  reflectance.  The  real  and 
imaginary  parts  of  the  Fourier  transform  of  the  object  centered  on 
the  optical  axis  are  shown  in  (A).  Since  the  object  is  symmetric  in 
this  case,  the  imaginary  part  of  the  transform  vanishes.  (B)  The 
object  was  translated  by  one- half  of  a  grating  cycle. 


transform  would  have  to  be  rotated  at  30  Hz,  corre¬ 
sponding  to  an  easily  obtainable  prism  rotation  rate  of 
7.5  Hz  =  450  rpm.  Indeed  much  higher  rates  have 
been  reported  with  excellent  image  quality.26 

Complex  transforms  obtained  with  this  system  at  a 
rate  of  2.5  frames/sec  are  presented  in  Figs.  3-6.  In 
each  case,  the  positive  part  of  the  Fourier  transform  is 
presented;  i.e.,  areas  of  the  transform  with  amplitude 
greater  than  zero  are  bright,  while  those  areas  with 
amplitude  less  than  zero  are  dark.  Note  the  difference 
in  the  usual  presentation  of  the  squared  magnitude  of 
the  Fourier  transform,  where  areas  with  amplitude 
both  greater  than  or  less  than  zero  are  bright,  and  the 
zero-crossings  are  dark.  In  Fig.  3,  the  object  was  a 
circular  aperture  1.5  mm  in  diameter.  In  the  first  case, 
the  aperture  was  centered  in  the  flying-line  scanner 
resulting  in  a  symmetric  object.  The  Fourier  trans¬ 
form  is  purely  real,  and  the  cosine  transform  is  the 
well-known  Airy  pattern.  If  the  aperture  is  translated 
in  the  scanner  so  that  the  object  is  no  longer  symmet¬ 
ric,  a  linear  phase  term  in  the  transform  appears  as 
fringes  in  the  output:  the  greater  the  shift,  the  larger 
the  frequency  of  the  linear  phase.  In  Fig.  4,  the  object 
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was  a  circular  aperture  over  a  transparent  grating 
whose  Fourier  transform  is  the  Airy  pattern  of  the 
circular  aperture  on  the  ±1  orders  of  the  grating. 
Each  transform  was  obtained  for  a  different  shift  of  the 
grating  in  the  aperture  resulting  in  different  phasing  of 
the  Airy  patterns.  In  Fig.  5,  the  object  was  a  pair  of 
circular  apertures,  each  1  mm  in  diameter  with  the 
centers  separated  by  5  mm.  The  first  case  is  the  cosine 
transform  of  the  centered  pair.  Note  the  phase  shift  in 
the  fringes  on  different  rings  of  the  Airy  pattern.  The 
other  cases  are  cosine  and  sine  transforms  of  the  pair 
with  the  optical  axis  positioned  1  mm  above  the  center 
line  through  the  apertures  and  with  the  axis  located  1 
mm  along  the  center  line. 

Figure  6  demonstrates  the  capability  of  this  system 
to  compute  the  complex  Fourier  transform  in  reflec¬ 
tion.  A  beam  splitter  was  inserted  into  the  system 
ahead  of  the  object  to  direct  the  reflection  of  the  flying 
line  onto  the  photomultiplier.  The  object  was  identi¬ 
cal  to  that  in  Fig.  4.  Since  the  overall  reflectance  is  less 
than  the  transmittance,  and  since  the  modulation  in 
reflectance  is  less  as  well,  the  SNR  of  the  reflective 
transforms  is  lower  than  that  in  the  transmissive  case. 
However,  the  capability  of  performing  complex  reflec¬ 
tive  transforms  with  this  system  is  clearly  demonstrat¬ 
ed. 

IV.  Conclusions 

We  have  demonstrated  a  hybrid  system  to  produce 
the  complex  2-D  Fourier  transform  in  nearly  real  time. 
The  transform  is  presented  as  real  and  imaginary  parts 
(i.e.,  cosine  and  sine  transforms)  in  separate  output 
channels.  Transmissive  or  reflective  functions  can  be 
transformed  in  this  manner. 
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Abstract 

It  is  well  known  that  many  mathematical  opera¬ 
tions  on  data  sets  of  dimension  two  or  higher  nay 
be  performed  by  reducing  the  data  to  one-dimen¬ 
sional  projections  by  means  of  the  Radon  transform. 
This  is  the  governing  principle  of  medical  computed 
tomography.  In  this  paper,  we  describe  a  system 
that  performs  the  Radon  transform  of  two-dimen¬ 
sional  images  and  uses  SAW  devices  to  perform  the 
data  processing.  Two  processing  operations  are  dem¬ 
onstrated:  Fourier  transf  ortna  tion  of  the  data  by 
means  of  the  chirp  transform,  and  convolution  of  the 
data  with  a  scored  filter  function  by  means  of  a 
SAW  correlator.  After  processing,  a  custom  SAW 
filter  and  an  optical  system  are  used  to  recon¬ 
struct  the  processed  image  in  two  dimensions.  The 
resolution  of  the  processor  is  currently  limited  by 
the  SAW  devices  (50  points  for  the  chirp  trans¬ 
former,  300  for  the  convolver),  but  better  devices 
are  available.  This  system  is  capable  of  performing 
two-dimensional  Fourier  transforms  at  video  rates 
(30  frames/s),  which  is  much  faster  than  current 
digital  systems.  An  extension  of  the  system  to  pro¬ 
cess  three-dimensional  data  is  described. 

In  trodmc  tion 

The  Radon  transform  has  received  much  atten¬ 
tion  in  the  scientific  community  since  the  invention 
of  x-ray  computed  tomography  (CT)  in  the  1960's.  It 
nas  found  application  in  such  diverse  disciplines  as 
astronomy,  nuclear  magnetic  resonance,  and  geophy¬ 
sics.  The  mathematics  of  the  transform  were  der¬ 
ived  and  published  by  Johann  Radon  in  1917  (i), 

where  he  proved  that  a  mathematical  function  can  be 
reconstructed  from  the  complete  set  of  its  line-in¬ 
tegral  projections.  In  the  case  of  CT,  measured  x- 
ray  transmissions  are  simply  related  to  the  line  in¬ 
tegral  of  the  x-ray  absorption  coefficient.  By 
taking  an  adequately  sampled  set  of  one-dimensiona  l 
data,  a  two-dimensional  map  of  the  x-ray  absorption 
coefficient  can  be  reconstructed,  usually  by  digital 
means. 

We  propose  to  use  the  Radon  transform  from  a 
different  perspective.  Instead  of  having  one-dimen¬ 
sional  projections  inherent  in  the  data  collection, 
we  use  the  Radon  transform  to  make  two-dimensional 
lata  susceptible  to  processing  by  fast  one-dimen¬ 
sional  devices.  Several  types  of  one-dimensional 
processors  exist;  the  SAW  filter  is  but  one.  Many 
two-dimensional  operations  can  be  performed  by 


means  of  the  Radon  transform,  e.g.,  spectral  ana¬ 
lysis,  convolution,  and  Fourier  filtering.  Such 
operations  can  be  performed  digitally,  of  course, 
but  the  process  may  be  time-consuming  and  the  pro¬ 
cessor  expensive.  By  operating  on  the  one-dimen- 
slonal  projection  Instead,  it  is  possible  that  the 
processor  may  be  significantly  faster  and/or  cheaper 
than  its  digital  counterpart.  Consider  Fourier 
filtering  of  a  two-dimensional  image,  for  example. 
Three  steps  are  required:  Fourier  transformation, 
filter  multiplication,  and  inverse  transformation. 
This  operation  may  be  performed  digitally,  by  coher¬ 
ent  optics,  or  with  one-dimensiona l  SAW  filters  by 
means  of  the  Radon  transform.  The  invention  of  the 
fast  Fourier  transform  (FFT)  algorithm  and  the  array 
processor  have  dramatically  speeded  up  digital  Four¬ 
ier  transform  calculations,  but  the  process  is  still 
slow.  A  typical  stand-alone  minicomputer,  the  DEC 
11/34,  taquiras  approximately  10  minutes  to  Fourier 
transform  a  512  x  512  9-bit  array.  Adding  an  array 
processor  speeds  this  by  an  order  of  magnitude  at 
significantly  increased  cost.  The  Cray- 1,  one  of  the 
fastest  digital  computers  ever,  still  requires  about 
l  second  to  perform  a  two-dimensional  Fourier 
transform  and  is  very  expensive.  Coherent  optics 
can  perform  Fourier  transforms  easily,  cheaply,  and 
at  the  speed  of  light,  but  the  output  Ls  noisy,  and 
there  are  still  no  fully  satisfactory  spatial  light 
modulators  to  allow  analysis  of  rapidly  time-vary¬ 
ing  inputs.  We  propose  to  perform  two-dimensional 
Fourier  transforms  by  operating  on  one-dimensional 
projections  with  SAW  chirp  filters.  The  resulting 
processor  should  be  inexpensive  relative  to  the  dig¬ 
ital  system,  but  more  importantly,  it  should  be 
fast:  we  envision  operation  at  video  rates. 

Thmory 

Mathematical  analysis  of  the  Radon  transform 
is  straightforward  and  has  been  treated  in  several 
references  (2,3).’  We  shall  touch  briefly  on  the  main 
points  relevant  to  the  application  at  hand,  ;.e., 
two-dimensional  Fourier  analysis  and  filtering. 

A  one-dimensional  projection  af  a  two- 

dimensional  function  f(r;  taken  aLong  azimuth  direc¬ 
tion  >  (relative  to  the  x  axis;  is  ietined  as 


d 2  r  f  (  r  ■  6  ( p  - 


r*n). 
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The  one-dimensiona  1  delta  function  reduces  the  area 
integral  to  a  Line  integral  along  a  line  at  an  angle 
*  to  the  x  axis  and  at  a  distance  p  from  the  origin 
(Figure  l).  The  set  ( X  ( ♦,  p )  1  for  all  azimuth  angles  * 
is  the  Radon  transform  of  f(r). 

By  caking  the  one-dioensiona  1  Fourier  transform 
of  a  line-integral  projection*  an  important  result  is 
obtained: 


F,[  X  ( »>p)l 


dp  exp(-2rlvp) 


d*r  f(r)d(p  -  p‘n) 


■  I  !  d'r  f(r)  exp(-2»invr)  -  F(p)|  , 
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where  capital  letters  denote  Fourier  transforms  of 
the  corresponding  lower-case  functions.  This  is  the 
central-slice  theorem.  In  words,  the  one-dimen¬ 
sional  Fourier  transform  of  a  projection  l*(p) 
yields  one  line  through  the  origin  of  the  two-dimen¬ 
sional  Fourier  transform  of  the  original  function 
f (r)  (Figure  l). 


Figure  -  Seometry  of  the  Radon  transform.  (a) 
Oerivation  of  one  projection  MS, pi  hv  line  integrals 
along  azimuth  angle  S.  ib)  Central  slice  theorem: 
the  one-dimensional  Fourier  transform  of  a  line-in¬ 
tegral  projection  yields  one  Line  through  the  two- 
dimensional  Fourier  transform  of  the  original  two- 
dimensional  function. 


3y  similar,  though  more  involved  reasoning,  it 
can  be  shown  that  the  original  function  f(r)  may  be 
reconstructed  from  the  projection  data  by  means  of 
the  inverse  Radon  transform 


[|N|A(«,v)l  exp(2»ivr-n)  (3) 


Again  In  words,  the  original  function  f(r)  may  be 
reconstructed  from  the  projections  [X(*,p)|  by:  (1) 
talcing  the  one-dimenslonal  Fourier  transform  of 
X(S,p);  (2)  multiplying  by  the  one-dimenslonal  filter 
|  v  |  (v-flltering);  (3)  taking  the  Inverse  one-dimen¬ 
sional  Fourier  transform;  (4)  smearing  the  function 
beck  over  the  originel  projection  direction  (this 
creates  a  two-dimensional  function  from  the  one- 
dlmenslonel  function  and  is  called  beck  projection); 
and  (3)  Integrating  over  *  (summation). 

If  we  multiply  the  one-dimensional  Fourier 
transform  data  A (»,v)  by  another  filter  function  as 
well,  the  reconstructed  function  Is  a  Fourier-fil¬ 
tered  version  of  f(r). 

Other  expressions  (and  hence  other  procedures 
for  taking  the  Inverse  Radon  transform  exist  and  are 
given  In  Reference  3. 

In  addition,  it  can  be  shown  that  by  convolving 
line-integral  projections  from  two  two-dimensional 
images,  and  reconstructing  by  the  procedure  of  Eq. 
(3),  the  resulting  two-dimensional  Image  is  the  con¬ 
volution  of  the  two  Input  Images. 

experiment 

We  constructed  a  system  using  the  Radon  trans¬ 
form  to  perform  two-dimensional  spectral  analysis 
using  SAW  chirp  filters.  The  apparatus  is  dia¬ 
grammed  in  Figure  3.  The  Radon  transform  of  the 
Input  transparency  is  derived  by  scanning  It  with  a 
line  of  HeNa  laser  light.  The  light  transmitted 
through  the  transparency  is  collected  on  a  photo¬ 
multiplier  tube  (PMT).  At  one  Instant  of  time,  the 
output  of  the  PUT  Is  proportional  to  the  line  integ¬ 
ral  of  the  Intensity  transmission  of  the  transpar¬ 
ency  along  the  line  of  light.  By  scanning  the  line 
perpendicular  to  Itself,  the  time  signal  from  the 
PMT  is  proportional  to  the  line-integral  projection 
along  one  azimuth.  Rotating  the  direction  of  scan 
allows  derivation  of  the  complete  3et  of  line-inte¬ 
gral  projections--the  Radon  transform.  For  obvious 
reasons,  this  device  is  called  a  flying-line  scanner 
(FLS). 

Recalling  the  centra  1 -s  lice  theorem,  we  know 
that  the  one-dimenslonal  Fourier  transform  of  one 
projection  yields  one  line  through  the  two-dimen¬ 
sional  transform  of  the  original  function.  By  using 
SAW  filters  in  the  chirp  Fourier  transform  algorithm 
(4,3,6),  the  Fourier  transform  of  each  projection  is 
taken  as  the  projection  data  are  derived.  We  used 
SAW  dispersive  filters  for  the  chirps  (Andersen  Labs 
models  OS- 120- L0-20-251A  and  -252A),  The  time 
dispersion  of  each  is  20  js  and  the  bandwidth  is  LO 
hHz.  The  time-bandwidth  product  T3W)  of  the  entire 
system  is  only  50,  but  filters  exist  that  could 
boost  this  to  2000  or  more. 


In  this  demonstration,  only  the  modulus  of  the 
Fourier  transform  is  computed,  but  we  plan  to  util¬ 
ize  a  third  chirp  filter  to  perform  the  post-nuiti- 
plication  and  derive  the  phase  information. 
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Figure  2  -  Diagram  of  a  two-dimensional  spectrum 

analyzer  using  SAW  devices.  The  projections  (X(*,p)] 
are  derived  by  the  flying- line  scanner.  The  one- 
dimensionai  Fourier  transforms  of  the  projections 
are  produced  by  the  SAW  chirp  filters.  The  Fourier 
transform  signal  modulates  the  CRT  trace.  The 
proper  azimuth  for  display  is  selected  by  the  image 
rota  cor. 


Extending  this  system  to  allow  complete  Four¬ 
ier  filtering  is  straightforward  and  is  diagrammed 
in  Figure  4.  The  Fourier  transforms  of  the  projec¬ 
tion  data  are  multiplied  by  a  filter  function,  which 
can  be  clocked  out  of  ROM  or  produced  by  a  function 
generator.  The  filtered  transforms  are  then  applied 
to  an  inverse  Radon  transformer  using  the  procedure 
of  Eq.  (3).  The  j  v  | -filter  multiplication  is  to  be 
done  using  a  custom  SAW  filter  that  is  presently 
being  constructed  in  the  University  of  Arizona 
Microelectronics  Laboratory.  The  inverse  one-dimen¬ 
sional  Fourier  transf orma tion  will  be  done  by  means 
of  the  chirp- transform  algorithm  and  the  output 
applied  to  the  z  axis  of  a  CRT.  To  perform  the 
back-projection  (smearing),  a  cylindrical  lens  is 
used  to  collimate  the  image  in  one  dimension.  The 
integration  over  azimuth  angle  is  carried  out  as 
before  by  rotating  the  image  on  the  recording  film. 
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To  complete  the  two-dimensional  spectral  ana¬ 
lysis,  it  is  necessary  to  display  the  transforms  of 
the  projections  at  the  proper  orientations.  After 
detection  and  amplification,  the  transform  of  the 
projection  is  applied  to  the  z  axis  of  a  CRT  whose 
trace  is  imaged  on  a  photographic  film.  As  the 
azimuth  of  scan  of  the  FLS  is  rotated,  the  image  of 
the  CRT  trace  is  rotated  at  the  same  rate,  building 
up  the  two-dimensional  Fourier  spectrum  modulus  on 
the  film.  A  result  from  this  experiment  is  shown  in 
Figure  3.  The  input  transparency  consisted  of  three 
gratings  oriented  at  various  angles:  two  fine 
crossed  gratings  overlaid  with  a-  section  of  coarse 
grating.  In  the  Fourier  transform  built  up  from  the 
projection  data,  the  fundamental  frequency  of  the 
fine  gratings  and  several  orders  of  the  coarse  grat¬ 
ing  are  visible.  This  spectrum  was  built  up  slowly, 
but  by  rapid  rotation  of  the  scan  direction,  we 
expect  to  perform  two-dimensional  spectral  analysis 
at  video  races  (30  frames/s)  or  faster. 


Figure  3  -  Results  of  two-dimensional  spectral  ana¬ 
lysis  using  SAW  devices. 

a)  Input  transparency  consisting  of  three  gratings. 
■  b;  Two-dimensional  spectrum,  showing  the  fundamen¬ 
tal  order  of  the  fine  gratings  and  several  orders  of 
the  coarser  grating. 


Figure  4  -  Slock  diagram  of  a  two-dimensiona  L 
system  to  do  Fourier  filtering.  The  signal  is  fil¬ 
tered  in  the  frequency  domain  and  transformed  beck 
to  the  space  domain  by  the  SAW  chirp- transf  orm 
algorithm.  The  cylindrical  lens  performs  the  beck- 
projection  (creates  a  two-dimensional  function  "it 
of  a  one-dimensional  function),  and  the  proper  azim¬ 
uth  is  selected  by  the  image  rotator. 


As  mentioned  previously,  another  operation  sus¬ 
ceptible  to  Radon  transform  analysis  is  the  convolu¬ 
tion  of  two  two-dimensional  functions.  The  neces¬ 
sary  apparatus  is  diagrammed  in  Figure  5.  Both 
inputs  may  be  projection  data  from  flying-line 
scanners,  but  it  is  often  useful  to  convolve  a  two- 
dimensional  function  with  a  scored  filter  function. 
This  function  may  be  scored  in  ROM,  clocked  out  to  i 
fast  D/A  converter,  and  used  to  modulate  a  carrier. 
The  resulting  signal  is  applied  to  one  input  of  3 
SAW  convolver.  The  projection  data  from  the  FLS 
modulate  the  carrier  and  are  applied  to  the  other 
input  of  the  convolver.  The  filter  function  may  be 
varied  with  azimuth  angle  by  clocking  a  different 
function  out  of  ROM  for  each  azimuth. 


Reconstruction  of  the  two-dimensional  convolu¬ 
tion  also  follows  the  procedure  of  Eq.  v3).  The 
function  is  v  I -filtered  in  the  custom  SAW  device, 
demodulated,  and  biased  jp  to  allow  display  at  bipo¬ 
lar  output.  This  signal  modulates  the  IRT  and  :s 
back-projected  and  integrated  over  the  azimuth  as 
before. 
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Figure  5  -  Block  diagram  of  a  system  to  perform 
two-dimensional  convolution  with  a  SaW  convolver. 
The  convolution  may  be  performed  between  two  two- 
dimensional  Inputs  or  between  a  two-dimensional 
input  and  a  stored  filter  function.  The  SAW  con¬ 
volver  output  is  jv|-filtered  by  a  custom  SAW 
filter. 
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Figure  6  -  Block  diagram  of  a  three-dimensiona 1 

spectrum  analyzer. 


over  the  digital  computer  example  above.  The  appli¬ 
cation  of  three-dimensional  processing  is  discussed 
further  in  Reference  7. 


Extension  to  Three-Dimensional  Date 

The  time-consuming  nature  of  three-dimensional 
data  processing  is  even  more  extreme  than  for  two- 
dimensional  data.  Performing  a  (512)*  FFT  on  a  min¬ 
icomputer  with  array  processor  and  fast  disk  memory 
may  take  two  days  or  more.  By  applying  the  princi¬ 
ples  of  the  Radon  transform,  we  expect  to  speed  up 
the  computation  considerably. 

In  the  three-dimensional  case,  the  Radon  trans¬ 
form  consists  of  the  complete  set  of  one-dimen¬ 
sional  integrals  taken  over  planes  of  the  three- 
dimensional  function.  The  three-dimensional  cen¬ 
tral-slice  theorem  states  that  the  one-dimensional 
Fourier  transform  of  a  planar  projection  yields  one 
Line  through  the  three-dimensional  Fourier  transform 
of  the  three-dimensional  function  (3). 

As  an  example,  consider  three-dimensional  spec¬ 
tral  analysis  of  a  function  stored  as  frames  of  a 
movie  film  (512  images,  each  5L2  x  512  pixels,  say). 
We  tan  use  SAW  chirp  filters  to  compute  the  three- 
dimensional  Fourier  transform.  A  block  diagram  is 
3hovn  *n  Figure  6.  The  data  manipulation  is  consid¬ 
erably  more  complicated  than  the  two-dimensional 
case,  since  the  Radon  transform  projections  are  now 
para  me  tenzed  by  two  angles.  But  by  using  a  digital 
video  frame  store  and  a  flying- Line  scanner,  we  can 
build  up  the  entire  Radon  transform,  sampled  at  512 
azimuth  angLes,  with  512  passes  of  the  movie  film. 
The  digital  frame  store  Is  then  read  out  through  a 
fast  D/A  to  the  5 AW  chirp- transformer.  The  projec¬ 
tion  transforms  modulate  the  CRT  as  before,  and  are 
imaged  onto  film.  We  build  up  the  512  frames  of  the 
transform  one  at  a  time  by  selecting  only  that  part 
of  the  projection  transform  relevant  to  the  frame 
at  hand.  After  reading  out  the  video  frame  store 
5 12  times,  the  compLete  three-dimensional  transform 
is  built  ip.  With  present  ’video  storage  technology, 
the  operation  is  envisioned  to  take  17  seconds  per 
frame,  or  Less  than  4  hours  for  the  complete  set. 
This  is  an  imorovement  of  an  order  of  magnitude 


Conclusions 

We  have  demonstrated  the  ability  of  one-dimen¬ 
sional  processing  devices,  such  as  SAW  filters,  to 
perform  certain  two-dimensional  processing  opera¬ 
tions  by  means  of  the  Radon  transform.  It  is  anti¬ 
cipated  that  this  will  allow  these  operations  to  be 
performed  much  more  rapidly  than  is  now  possible 
with  digital  techniques. 
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Abstract.  The  well-known  central-slice,  or  projection-slice,  theorem  states  that 
the  Radon  transform  can  be  used  to  reduce  a  two-dimensional  Fourier  trans¬ 
form  to  a  series  of  one-dimensional  Fourier  transforms.  In  this  paper  we 
describe  a  practical  system  for  implementing  this  theorem.  The  Radon  trans¬ 
form  is  carried  out  with  a  rotating  prism  and  a  flying-line  scanner,  while  the 
one-dimensional  Fourier  transforms  are  performed  with  surface  acoustic  wave 
filters.  Both  real  and  imaginary  pans  of  the  complex  Fourier  transform  can  be 
obtained.  A  method  of  displaying  the  two-dimensional  Fourier  transforms  is 
described,  and  representative  transforms  are  shown.  Application  of  this 
approach  to  Labeyrie  speckle  interferometry  is  demonstrated. 

Subject  terms:  optical  computing:  Radon  transform;  Fourier  transform;  speckle 
interferometry. 
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1.  INTRODUCTION 

The  two-dimensional  Radon  transform  reduces  a  2-D  function  to  a 
senes  of  l-D  functions  by  integrating  over  a  senes  of  lines.  Although 
this  tran'^nn  is  best  known  in  connection  with  image  reconstruc¬ 
tion  from  projections,  as  in  medical  computed  tomography,  it  is  also 
useful  in  general  signal-processing  or  image-processing  applications. 
Many  operations  that  can  be  performed  on  a  2-D  function  can  also 
be  done  by  performing  l-D  operations  on  the  projections.  Recent 
work  has  demonstrated  the  usefulness  of  this  approach  in  calculating 
Fourier  transforms 1  ;  and  Wigner  distribution  functions,'  as  well  as 
in  pattern  recognition,4  s  image  filtering,*  7  and  bandwidth 
compression.4 

That  these  operations  are  possible  in  the  l-D  Radon  domain  is  a 
consequence  of  the  celebrated  central-slice,  or  projection-slice, 
theorem.  This  theorem  states  that  if  a  l-D  projection  of  a  2-D 
function  is  formed  by  integrating  over  a  set  of  parallel  lines,  the  !-D 
Founer  transform  of  the  projection  is  one  line  through  the  2-D 
Fourier  transform  of  the  function  itself  (see  Fig.  1 ).  This  line  passes 
through  the  origin  of  the  2-D  Founer  space  ( hence  the  term  central 
slice).  By  varying  the  orientation  of  the  lines  of  integration,  the  whole 
2-D  Founer  space  can  be  mapped  out  in  a  polar  format. 

In  this  paper  we  describe  in  detail  a  practical  system  for  perform¬ 
ing  2-D  Founer  transforms  in  the  Radon  domain.  Special  attention 
is  given  to  the  electronics  for  displaying  the  2-D  Founer  transforms, 
and  several  representative  transforms  are  shown.  As  an  illustration 
of  this  approach,  we  demonstrate  that  the  Radon  transform  can  be 
used  to  process  data  from  astronomical  speckle  interferometry. 

Invited  Piper  OP- 1 1 1  received  June  1 1,  >984:  accepted  for  publication  July  30.  1984. 
received  by  Managing  Editor  Sept.  24.  1984 
•  1985  Society  ot  Photo-Optical  Initrumemation  Engineers. 


2.  PRODUCING  AND  TRANSFORMING  THE  ONE¬ 
DIMENSIONAL  DATA 

In  this  section  we  describe  the  subsystems  for  producing  the  Radon 
transform  and  the  l-D  Fourier  transforms.  Since  these  subsystems 
have  been  described  previously,1  -J  only  a  brief  review  is  given  here. 

Assume  that  the  2-D  function  to  be  transformed  is  in  the  form  of  a 
photographic  transparency  or  print.  The  projection  data  A^lp)  are 
derived  from  the  input  function  f(x,y)  by  scanning  a  line  of  light 
perpendicular  to  itself  across  the  function  at  an  angle  <f>  (see  Fig.  1 ). 
The  perpendicular  distance  of  the  line  from  the  origin  is  p.  and  the 
line  is  uniquely  specified  by  the  variables  p  and  <b.  The  light  transmit¬ 
ted  or  reflected  by  the  object  is  detected  by  a  photomultiplier  tube 
(PMT).  The  signal  out  of  the  PMT  is  then  proportional  to  the  line 
integral  of  the  object  transmittance  or  reflectance  along  the  line 
(p ,  ib).  As  the  line  is  scanned  by  means  of  an  acousto-optic  deflector, 
the  variable  p  changes,  and  one  scan  produces  one  projection  p). 
A  rotating  prism  in  the  system  changes  the  orientation  of  the  line, 
which  is  always  scanned  perpendicular  to  itself,  and  provides  other 
projections  in  the  data  set.  In  this  way  the  entire  data  set.  sampled  in 
<t>  but  continuous  in  p,  can  be  formed.  This  system  is  referred  to  as  a 
flying-line  scanner.1-' 

The  l-D  Fourier  transforms  A^lv)  are  formed  by  a  surface 
acoustic  wave  (SAW)  chirp  transformer'-9  in  which  the  input  signal 
(the  projection)  is  premultiplied  by  a  chirp  produced  by  impulsing  a 
SAW  device.  The  resulting  signal  is  filtered  (convolved)  by  a  second 
SAW  chirp  filter  in  which  the  chirp  rate  is  equal  and  opposite  to  that 
of  the  premultiply  signal.  The  signal  out  of  this  second  filter  is 
coherently  detected  with  a  third  chirp  as  a  reference.  In-phase  and 
quadrature  outputs  of  the  coherent  detector  give,  respectively,  the 
real  and  imaginary  parts  of  the  complex  Fourier  transform.  If  only 
the  modulus  of  the  transform  is  desired,  the  third  chirp  can  be 
omitted  and  incoherent  detection  used. 

By  the  central-slice  theorem.  A^( p)  is  also  the  2-D  Fourier  trans¬ 
form  of  f(x.y)  evaluated  at  polar  coordinates  (p,d>)  in  the  2-D 
Fourier  space,  where  p  =  |  v  | . 

3.  DISPLAYING  THE  TRANSFORM 

A  simple  way  to  display  the  l-D  Fourier  transforms  is  in  the  so-called 
“sinogram"  format,  in  which  the  radial  frequency  variable  p  is  plot¬ 
ted  horizontally  and  the  azimuthal  variable  <b  is  plotted  vertically. 
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This  representation  is  certainly  legitimate  and  useful  in  some  applica¬ 
tions.*  but  it  is  often  desirable  to  present  the  data  in  polar  format. 
This  not  only  makes  the  transform  more  recognizable  to  someone 
not  familiar  with  sinograms  but  also  presents  the  data  in  a  form  that 
can  be  further  processed  in  cascaded  systems.  A  system  to  accom¬ 
plish  this  polar  display  has  been  designed,  built,  and  operated  and  is 
described  in  this  section. 

For  now.  assume  that  only  the  modulus  of  the  2-D  transform  is  to 
be  displayed.  The  rf  signal  from  the  SAW  chirp  transformer  is 
detected  incoherently,  forming  a  signal  proportional  to  the  squared 
magnitude  of  the  Fourier  transform,  which  is  then  used  to  intensity- 


Fig.  1.  (a)  Geometry  of  the  Radon  tranaform  and  (b|  illustration  of  the 
central-slice  theorem. 


modulate  a  spot  on  a  CRT  display  that  is  being  scanned  in  a  polar 
raster!  Fig.  2).  The  polar  angle  of  the  raster  is  the  same  as  the  angle  o 
that  specifies  the  orientation  of  the  line  of  light  in  the  living-line 
scanner,  while  the  radial  variable  on  the  raster  corresponds  to  the 
frequency  p.  The  time-averaged  intensity  on  the  screen  represents  the 
2-D  Fourier  transform  in  a  direct  format  and  is  equivalent  to  the 
intensity  distribution  in  the  Fourier  plane  of  a  coherent  optical 
transformer  with  the  same  input  function. 

To  maintain  synchronism  between  the  polar  raster  and  the  living- 
line  scanner,  a  stepper  motor  is  used  to  control  the  rotation  of  the 
prism  in  the  scanner.  Each  step  of  the  stepper  motor  changes  the 
orientation  of  the  scanning  line  by  tr  100  rad.  A  free-running  circuit 
operates  the  stepper  motor  from  about  1 , 2  to  1 000  steps  per  second . 
Each  time  a  step  occurs,  a  short  transistor-transistor  logic  (TTL) 
pulse  is  sent  to  a  Commodore  64  computer.  Upon  receiving  this 
pulse,  the  computer  updates  an  index  register  to  indicate  the  new 
angle,  sends  bytes  representing  the  sine  and  cosine  of  the  new  angle  to 
two  digital-to-anaiog  converters,  and  finally  sends  a  short  TTL  pulse 
to  a  third  output  port  to  signal  the  rest  of  the  system  to  generate  a  new 
line  of  data. 

The  start-of-line  pulse  from  the  computer  starts  a  scan  in  the 
flying-line  scanner,  triggers  the  impulse  generator  for  the  premultiply 
chirp,  and  triggers  a  delay  circuit  whose  output  after  the  proper  delay 
is  a  30  ps  pulse  used  to  control  the  display  The  delay  is  adj usted  such 
that  the  Fourier  transform  data  are  centered  within  the  30  ms  pulse. 
During  this  pulse,  a  bipolar  ramp  function  is  generated,  passing 
through  zero  at  the  same  time  the  zero-frequency  component  of  the 
Fourier  transform  is  available.  This  ramp  function  is  multiplied  by 
the  sine  and  cosine  values,  and  the  results  of  these  multiplications  are 
used  to  control  the  x  and  y  deflections  of  a  spot  on  a  CRT  display 
This  causes  the  spot  to  travel  across  the  screen  at  a  constant  speed  at 
an  angle  equal  to  the  scan  angle  <p.  reaching  the  center  of  the  screen  at 
a  time  corresponding  to  the  zero-frequency  output  time  of  the 
Fourier  transformer. 

If  the  signal  coming  out  of  the  SAW  transformer  were  simply 
detected  and  used  to  intensity-modulate  the  CRT.  the  screen  would 
display  the  desired  output  except  for  one  problem.  As  the  entire 
output  is  built  up.  the  radial  scanning  pattern  fills  the  space  near  the 
center  much  more  densely  than  near  the  edges.  The  resulting  time- 
averaged  intensty  distribution  would  appear  as  Flp)  times  I  o. 
where  p  is  the  2-D  frequency  vector,  p  is  its  magnitude,  and  Fl  p)  is 
the  2-D  Fourier  transform  of  the  input  function  fix.  v).  I  Recall  that  p 
is  also  the  magnitude  of  the  l-D  frequency  v.  but  v  can  be  bipolar, 
while  p  is  always  nonnegative.)  In  order  to  eliminate  the  I  p  weight¬ 
ing,  it  is  necessary  to  multiply  the  signal  before  detection  by  \  v  This 
is  accomplished  in  the  following  manner  The  ramp  function  driving 
the  multipliers  is  used  as  the  input  to  an  absolute-salue  amplifier 
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Fig.  3.  2-0  Fourier  transform  of  a  square- wave  grating  in  a  small  circular 
aperture. 


Fig.  4.  2-0  Fourier  transform  of  two  crossed  gratings  making  an  angle  of 
approximately  90° . 


The  output  of  this  amplifier  is  a  V-shaped  function  that  is  the  desired 
multiplier. 

Because  the  signal-to-noise  ratio  is  generally  decreasing  with 
increasing  frequency,  it  is  also  desirable  to  roll  off  or  apodize  the  V 
function  at  higher  frequencies.  This  is  easily  accomplished  by  current- 
limiting  the  output  of  the  V  generator.  The  signal  resulting  from 
multiplying  the  output  of  the  convolving  filter  by  the  apodized  V 
function  is  square-law  detected  and  used  to  intensity-modulate  the 
spot  on  the  CRT.  The  resulting  display  from  one  scan  is  a  line  in  2-D 
Fourier  space,  filtered  by  pand  the  apodizing  function.  As  all  angles 
are  traced  out.  an  entire  disk  of  Fourier  space  is  built  up  on  the 
screen. 

It  is  straightforward  to  extend  this  system  to  a  CRT  display  of 
complex  Fourier  transforms.  The  coherent  detector  in  the  chirp 
transformer  provides  bipolar  signals  proportional  to  the  real  and 
imaginary  parts  of  the  complex  transform.  These  signals  can  be 
separated  further  into  four  nonnegative  signals,  namely,  the  positive- 
real.  negative-real,  positive-imaginary,  and  negative-imaginary 
components,  each  of  which  can  be  used  to  intensity-modulate  a  CRT 
display.  Either  four  separate  CRTs  can  be  used,  or  a  single  display 
can  be  used  sequentially  for  the  four  components.  Alternatively, 
analog  electronic  modules  are  available  to  convert  the  real  and 
imaginary  parts  to  modulus  and  phase,  which  can  be  displayed  with 
the  system  described  above. 

4.  RESULTS 

Several  examples  of  2-D  Fourier  transforms  produced  on  the 
Radon-Founer  transformer  are  shown  in  Figs.  3  through  7. 

Figure  3  shows  the  transform  of  a  square-wave  grating  with  a  duty 
cycle  of  about  0.7;  the  ±  1  and  ±2  orders  are  seen.  The  aperture  of  the 
grating  is  a  small  circular  iris,  and  the  rings  of  the  Airy  disk  are  visible 
in  the  ±  1  orders.  The  effects  of  the  angular  sampling  can  be  seen  in 
the  ±2  orders  since  only  four  or  five  sweeps  of  the  flying-line  scanner 
intersect  these  orders. 

Figure  4  shows  the  transform  of  two  overlapping  orthogonal 


Fig.  S.  2-0  Fourier  transtorm  of  two  croxaxd  grating*  making  an  angla  of 
approximately  48° . 


Fig.  6.  2-0  Fourier  transform  of  a  computer-generated  hologram  with  an 
annular  impulse  response. 


Fig.  7.  2-0  Fourier  transform  of  a  double  pinhole. 


square-wave  gratings.  Again,  the  duty  cycle  of  the  gratings  is  about 
0.7.  The  product  orders  resulting  from  the  convolution  ot  ihe  two 
individual  grating  spectra  are  clearly  seen.  Figure  5  is  similar  except 
that  the  angle  between  the  gratings  is  approximately  45° 

Figure  6  shows  the  transform  of  a  computer-generated  hologram 
that  has  an  annular  impulse  response. 

Figure  7  is  the  transform  of  a  double  pinhole  The  two  notches  are 
due  to  limited  dynamic  range  in  the  rf  mixers,  a  problem  that  can  be 
solved  with  better  mixers. 

With  the  SAW  filters  we  actually  used,  the  time-bandwidth  pro¬ 
duct  in  the  I -D  transforms  was  only  about  50.  so  the  results  shown  in 
these  figures  have  relatively  low  resolution,  containing  roughly  2000 
resolvable  spots.  However,  this  is  by  no  means  a  fundamental  limita¬ 
tion;  SAW  filters  are  commercially  available  that  will  provide  6000 
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Fig.  8.  SpacU*  intarfarogram  of  a  aimulatad  binary  (tar.  Tha  original 
speckle  pattama  ware  producad  on  a  computar.  but  thia  figure  was  peo- 
ducad  by  tha  Radon-Fouriar  transformer. 

spots  in  a  l-D  Fourier  transform,  or  almost  30  million  in  a  2-D 
transform. 

For  our  experiments,  the  time  required  to  produce  a  complete 
2-D  Fourier  transform  was  about  0.3  s.  but  again  this  is  not  a 
fundamental  limitation;  rather  it  is  a  limitation  on  how  rapidly  the 
prism  in  the  flying-line  scanner  could  be  rotated.  It  is  easily  possible 
to  rotate  a  prism  at  450  rpm,  which  would  yield  2-D  transforms  at 
video  rates,  and  even  I  ms  per  transform  appears  feasible. 

5.  SPECKLE  INTERFEROMETRY 

Astronomical  speckle  interferometry  is  an  ingenious  technique 
invented  by  Labeyrie  to  obtain  diffraction-limited  resolution  from  a 
telescope  in  spite  of  phase  perturbations  by  the  atmosphere. 10  In  this 
technique,  a  series  of  photographic  exposures  is  made,  with  each 
exposure  time  being  short  compared  to  the  scintillation  time  of  the 
atmosphere.  Each  image  is  Fourier  transformed,  either  digitally  or 
optically,  the  sum  of  the  squared  moduli  of  the  Fourier  transforms  is 
accumulated.  One  final  Fourier  transform  then  yields  the  autocorre¬ 
lation  of  the  object  with  diffraction-limited  resolution. 

Since  this  method  involves  a  large  number  of  Fourier  transforms, 
it  is  natural  to  consider  the  use  of  the  Radon  transform  to  reduce  the 
2-D  Fourier  transforms  to  l-D.  Indeed,  in  some  infrared  applica¬ 
tions.  I  -D  projections  of  speckle  patterns  are  observed  directly  by  use 
of  a  scanning  slit  in  the  image  plane."12 

To  demonstrate  the  use  of  our  Radon-Fourier  transformer  in 
speckle  interferometry,  we  simulated  a  series  of  20  speckle  patterns 
on  a  digital  computer.  Each  speckle  pattern  consisted  of  SO  pairs  of 
ellipses  of  random  size  and  ellipticity  but  with  constant  spacing 
between  members  of  the  same  pair.  The  resulting  configuration  of 
ellipses  was  intended  to  represent  the  speckle  pattern  that  would  be 
produced  by  a  binary  star.  The  20  speckle  patterns  were  photo¬ 
graphed  on  35  mm  transparency  film,  and  after  development  the  film 
strip  was  pulled  through  the  input  plane  of  the  Radon-Fourier 
transformer.  The  modulus  of  the  2-D  Founer  transform  was  dis¬ 
played  o.i  a  CRT  as  described  above,  and  a  camera  with  an  open 
shutter  was  used  to  accumulate  the  sum  of  the  Founer  moduli.  The 
resulting  image.  Fig.  8.  clearly  shows  the  fringe  pattern  characteristic 
of  a  double  star.  One  further  2-D  Founer  transform,  also  earned  out 
with  the  Radon-Fourier  transformer,  yielded  the  autocorrelation  of 
the  double  star,  as  shown  in  Fig.  9. 


Fig.  9.  Fourier  transform  of  Fig.  8.  which  it  tha  autocorrelation  function 
of  tha  simulated  binary  star. 


6.  CONCLUSIONS 

We  have  shown  that  the  Radon  transform  is  a  convenient  and  rapid 
vehicle  for  the  calculation  of  2-D  Fourier  transforms.  The  particular 
system  described  here,  which  is  based  on  a  flying-line  scanner  and  a 
SAW  chirp  Fourier  transformer,  has  a  number  of  advantages  over 
coherent  optical  Fourier  transformers.  It  does  not  require  that  the 
function  to  be  transformed  be  in  the  form  of  a  transparency;  it  works 
also  when  the  function  is  recorded  as  a  photographic  print  or  is  a 
natural  reflecting  scene.  Although  it  uses  a  laser  as  a  convenient 
source,  its  operation  does  not  depend  on  the  coherence  of  the  source. 
Furthermore,  the  full  complex  Fourier  transform  is  available,  some¬ 
thing  that  is  very  difficult  to  obtain  with  coherent  optical  techniques. 
The  system  is  also  extremely  fast.  With  presently  available  SAW 
filters,  a  system  similar  to  the  one  described  here  could  be  built  that 
would  produce  a  500  X  500  (500  points  across  the  diameter  of  the 
Fourier  plane  and  500  angles  in  the  range  0  to  n-)  2-D  transform  in 
1  /  30  s  and  a  5000  X  5000  transform  in  a  few  seconds. 
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A  bandwidth-compression  scheme  for  two-dimensional  data  is  presented  that  incorporates  the  Radon  transform. 
There  are  three  advantages  to  this  approach:  only  one-dimensional  operations  are  required,  the  dynamic  range 
requirements  of  the  compression  are  reduced  by  a  filtering  step  associated  with  the  inverse  Radon  transform,  and 
the  technique  is  readily  adaptive  to  the  data  structure.  A  rectilinear  object  is  compressed  to  demonstrate  the  algo¬ 
rithm. 


Introduction 

The  Radon  transform1*3  is  best  known  as  the  theoretical 
backbone  of  computed  tomography,  the  technique  that 
produces  cross-sectional  maps  of  x-ray  attenuation. 
This  transform  entails  projecting  a  two-dimensional 
slice  of  an  object’s  x-ray  attenuation  coefficient  along 
a  given  direction  in  the  plane  of  the  slice,  forming  a 
one-dimensional  data  set  for  each  projection  direction. 
Thus  the  two  spatial  dimensions  of  the  slice  are  trans¬ 
formed  into  one  spatial  and  one  angular  dimension  in 
Radon  space.  This  reduction  of  spatial  dimension  can 
be  used  to  reduce  two-dimensional  operations  on  a 
two-dimensional  object  to  a  set  of  one-dimensional 
operations  on  one-dimensional  objects,  with  each 
member  of  the  set  corresponding  to  a  projection  angle. 
In  particular,  as  a  consequence  of  the  central-slice 
theorem,  the  Radon  transform  makes  the  two-dimen¬ 
sional  Fourier  transform  of  a  two-dimensional  function 
readily  accessible  without  two-dimensional  operations’ 
actually  having  to  be  performed.  Thus  the  motivation 
exists  for  exploring  the  use  of  the  Radon  transform  in 
areas  outside  clinical  tomography.45  A  particularly 
direct  application  is  to  bandwidth  compression.6 

Compressing  the  data  necessary  to  represent  an 
image  (with  minimum  image  degradation)  is  important 
for  two  reasons:  storage  requirements  are  reduced,  and 
transmission  bandwidth  requirements  are  reduced.  If 
we  define  the  data  set  to  be  an  image  of  N  X  N  pixels 
with  each  pixel  corresponding  to  M  gray  levels,  com¬ 
pression  can  be  imposed  in  the  spatial  domain  or  in  a 
transform  domain.  Spatial  compression  consists  of 
reducing  (quantizing)  the  number  of  gray  levels  per 
pixel  and/or  reducing  the  number  of  pixels  in  the  image 
( i.e.,  reducing  the  radiometric  and  spatial  redundancy, 
respectively).  Transform  compression  consists  of 
transforming  the  image  (e.g.,  Fourier,  Hadamard,  Haar) 
and  then  quantizing  and/or  eliminating  the  coefficients 
of  the  transformed  image.7  To  reconstruct  the  image, 
the  inverse  transform  of  the  compressed  coefficients  is 
taken. 

The  Radon  transform  lends  itself  to  Fourier-trans- 
form  compression  for  three  reasons: 


(1)  The  entire  coding  process  can  be  performed  with 
state-of-the-art  one-dimensional  devices. 

(2)  The  large  dynamic  range  typical  of  the  compo¬ 
nents  of  the  Fourier  transform  is  significantly  reduced 
by  the  filtering  operation. 

(3)  One  line  through  the  center  of  the  two-dimen¬ 
sional  Fourier  transform  can  be  examined  at  a  time  and 
adaptively  compressed. 

Theory 

Radon  Transform 

The  Radon  transform  and  its  inverse2-3  are  central  to 
the  compression  technique.  Given  a  two-dimensional 
function  f(r),  where  r  is  the  spatial-position  vector  (x, 
y),  the  set  of  one-dimensional  projections  of  f  along  a 
given  direction  <t>  can  be  written  as 

X®(p)  =  /(r)<5(p  -  r  •  n)d2r.  (1) 

where  5(p  —  r  •  n)  is  a  one-dimensional  Dirac  delta 
function  restricting  the  integration  of  /  to  a  line  iwith 
normal  n)  located  a  distance  p  from  the  origin.  Thus, 
for  each  projection  direction  o.  a  one-dimensional 
function  X«(p)  is  constructed.  The  set  of  all  \,,lp  H — ®> 
<p<°°,0<tf><<r)  constitutes  the  Radon  transform 
of  f(x,y). 

Performing  the  one-dimensional  Fourier  transform 
on  Eq.  (1)  and  using  the  sifting  property  of  the  delta 
function  results  in 

A*(i')  =  jy_  /(r)  exp(-2xti'r  •  n)d-’r  =  Ftp) | 

(2) 

where  v  is  the  frequency-variable  conjugate  to  p.  p  is  the 
frequency- variable  conjugate  to  r.  and  Ftp)  is  the  two- 
dimensional  Fourier  transform  of  / ( r ),  evaluated  along 
the  line  p  =  ny,  where  n  in  the  frequency  domain  is 
parallel  to  n  in  the  spatial  domain.  This  is  the  cen¬ 
tral-slice  theorem. 

By  writing  fir)  in  terms  of  its  inverse  Fourier  trans- 


0146-9592/83/070395-03$  1 .00/0  ©  1983,  Optical  Society  of  America 


form  in  polar-coordinate  form  and  using  Eq.  (2),  we  can 
write  the  inverse  Radon  transform  as 

/( r)  =  f  dd>  [  P  d('j|t'|  exp(2irii'p) 

JO  [J p-n-r 

(3) 

If  one  looks  at  the  bracketed  term  in  Eq.  (3),  the  fol¬ 
lowing  operations  are  evident:  The  one-dimensional 
Fourier  transform  of  each  Radon  projection  A*(i>)  is 
multiplied  by  the  frequency  filter  |  x/| ;  then  the  inverse 
one-dimensional  Fourier  transform  is  applied  to  this 
product  and  evaluated  at  p  =  n  •  r.  This  is  the 
backprojection  step.  The  integral  over  <t>  is  the  sum¬ 
mation  of  the  backprojections  to  produce  f(r). 

Compression  Scheme 

From  the  central-3lice  theorem  of  Eq.  (2),  access  to  a  line 
passing  through  the  center  of  the  two-dimensional 
Fourier  transform  is  immediately  available  by  the 
one-dimensional  Fourier  transform  of  a  given  Radon 
projection.  This  suggests  that  an  adaptive  trans- 
form-compression  scheme  can  be  applied  to  one  line  of 
the  two-dimensional  Fourier  transform  at  a  time.  In 
fact,  the  compression  step  can  be  advantageously  ap¬ 
plied  to  the  filtered  line,  i.e.,  |v|  A0(t/)  of  Eq.  (3).  The 
filtered  and  compressed  line  is  then  stored  or  trans¬ 
mitted.  To  reconstruct  the  image,  the  inverse  Fourier 
transform  is  applied  to  each  previously  compressed  line, 
the  result  is  backprojected,  and  the  backprojections  are 
summed  to  produce  the  final  image. 

The  compression  is  accomplished  by  thresholding 
and  quantizing  the  components  of  each  Fourier  line. 
Because  the  projection  X*(p)  is  real,  its  Fourier  trans¬ 
form  A0(g)  is  hermitian  (i.e.,  the  real  part  is  even;  the 
imaginary  part  is  odd),  so  that  only  the  positive  half  {v 
>  0)  of  each  Fourier  line  need  be  transmitted  or  stored. 
The  thresholding  that  we  apply  is  to  truncate  each  line 
past  some  cutoff  frequency  C®,  which  is  variable  from 
line  to  line  (i.e.,  depends  on  <t>).  The  value  of  C,a  is  found 
from 

f  *  ||t<|  A*(p)  dv  =  S*  Sf-  )  •  (4) 

JO  |  |  Warn  ax/ 

where 

S0  =  |HA*(")  dp,  (5) 

S.tauu  is  the  largest  value  of  S*  for  0  <  <t>  <  x,  and  T  is 
a  parameter  that  controls  the  degree  of  truncation. 
Note  that  the  line  corresponding  to  S<,m u  is  never 
truncated  and  that,  as  T  -*  0,  C*  —  <*»  for  all  lines  (limit 
of  no  truncation).  This  method  of  choosing  C*  is  not, 
of  course,  fundamental:  other  algorithms  may  be  de¬ 
rived. 

After  truncating  the  line,  we  quantize  the  components 
by  dividing  the  full  dynamic  range  (positive  to  negative) 
specific  to  the  line  into  a  series  of  uniform,  discrete 
ranges.  Actually,  two  dynamic  ranges  exist,  one  each 
for  the  real  and  imaginary  parts.  The  component  that 
falls  within  a  particular  range  is  assigned  the  constant 
value  for  that  range. 


The  advantages  of  the  Radon  approach  are  now  dis¬ 
cussed. 

From  an  implementation  point  of  view,  hardware 
devices  for  carrying  out  one-dimensional  operations  are 
well  developed.  The  operations  for  each  projection  at 
the  compression  end  involve  a  one-dimensional  Fourier 
transform,  multiplication  by  a  linear  filter,  thresholding, 
quantizing,  and  coding  for  transmission  or  storage.  At 
the  receiving  end,  only  a  one-dimensional  inverse 
Fourier  transform  is  required,  followed  by  the 
backprojection  operation. 

The  dynamic  range  of  a  line  through  the  center  of  a 
two-dimensional  Fourier  transform  is  large.  To 
quantize  such  a  range  efficiently,  a  variable  quantizer 
would  be  required.  Multiplying  by  the  1 1/|  filter,  how¬ 
ever,  reduces  the  dynamic  range  of  the  line  (near  \i>\  = 
0,  where  the  components  are  usually  largest),  simpli¬ 
fying  the  requirements  of  the  quantizer. 

The  third  advantage  is  related  to  the  image-depen¬ 
dent  adaptability  of  the  compression  scheme.  An 
image  with  relatively  sharp,  straight  edges,  oriented  in 
particular  directions,  will  exhibit  a  transform  with  the 
energy  distributed  along  conjugate  directions,  de¬ 
pending  on  the  symmetry  of  the  original  image.  Be¬ 
cause  the  Radon  transform  handles  one  Fourier- 
transform  line  at  a  time,  each  filtered  line  can  be  ad¬ 
aptively  compressed  to  take  advantage  of  the  structure 
in  the  two-dimensional  Fourier-transform  plane. 


(c)  U) 


Fig.  1.  (a)  Reconstruction  from  Radon  projections  without 

thresholding  and  nominal  8-bit  quantization  (8  bits/pixel), 
(b)  truncation  of  48%  of  components  with  3-bit  quantization 
(1.6  bits/pixel),  (c)  truncation  of  66%  of  components  with  3-bit 
quantization  (1.1  bits/pixel),  (d)  truncation  of  48%  of  com¬ 
ponents  with  two-bit  quantization  (1.1  bits/pixei). 
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Fig.  2.  Truncated  and  quantized  Fourier  components  of  Fig. 
1(d). 


Illustrative  Example 

To  demonstrate  the  method,  the  object  illustrated  in 
Fig.  1(a)  is  compressed.  The  region  to  be  compressed 
is  a  circle  with  a  radius  of  64  pixels,  yielding  an  area  of 
12,868  pixels.  The  rectilinear  nature  of  the  figure  is 
useful  in  demonstrating  the  variable  compression  with 
projection  angle. 

The  object  is  viewed  by  a  TV  camera  through  an 
image  rotator.  The  one-dimensional  Radon  projection 
<\*(p)  of  the  object  is  obtained  (at  the  angle  <t>  defined 
by  the  image  rotator)  by  summing  the  camera  output 
along  a  horizontal  raster  line,  giving  the  digital  equiv¬ 
alent  of  the  line  integral  of  the  object  along  the  line. 
This  summation  is  performed  for  each  of  128  raster  lines 
to  yield  a  128-element  projection.  We  then  rotate  the 
image  rotator  through  a  small  angle  (1.8  deg)  and  find 
the  next  projection,  until  a  total  of  100  projections,  each 
with  128  elements,  has  been  taken. 

Each  projection  is  Fourier  transformed,  yielding  a  line 
through  the  center  of  the  two-dimensional  transform 
of  the  original  object.  This  line  is  then  filtered  [mul¬ 
tiplied  by  M;  see  Eq.  (3)]  and  compressed  by  thresh¬ 
olding  and  uniformly  quantizing  the  components,  as 
described  above.  The  object  is  reconstructed  by  taking 
the  one-dimensional  inverse  Fourier  transform  of  the 
filtered  and  compressed  line,  backprojecting,  and 
summing  over  all  projection  angles.  Figure  1(a)  illus¬ 
trates  the  reconstruction  from  the  Radon-transformed 
original  object  without  thresholding  or  quantizing  to 
provide  a  control  case. 

The  measure  of  the  compression  in  bits  per  pixel  (the 
bit  rate)  is  determined  by  counting  the  total  number  of 
bits  required  to  3tore  or  transmit  the  image  (including 
any  overhead)  divided  by  the  number  of  pixels  in  the 
image.  Figure  1(a)  has  approximately  8.0  bits/pixel 
because  there  are  100  angles  times  128  Fourier  compo¬ 
nents  per  angle,  times  8  bits  per  Fourier  component, 
divided  by  12,868  pixels. 

Three  different  compressions  of  Fig.  1(a)  are  now 


demonstrated.  Figure  1(b)  represents  thresholding 
with  T  =  0.3  and  quantization  of  the  Fourier  compo¬ 
nent’s  full  range  to  eight  gray  levels  (3  bits).  The 
thresholding  eliminates  48%  of  the  Fourier  components. 
The  overall  bit  rate  is  1.6  bits/pixel.  Figure  1(c)  has  T 
=  0.6,  eliminating  66%  of  the  components,  again  with 
eight  gray  levels  per  component,  giving  a  bit  rate  of  1.1 
bits/pixel.  To  investigate  a  coarser  quantization,  Fig. 
1(d)  represents  the  same  thresholding  as  in  Fig.  1(b)  but 
with  four  gray  levels  per  component,  giving  a  bit  rate  of 
1.1  bits/pixel.  Figure  2  is  a  representation  of  the 
truncated  and  quantized  components  that  produce  Fig. 
1(d).  The  adaptive  nature  of  the  compression  is  evident 
for  this  type  of  object.  Note  that  an  overhead  of  ap¬ 
proximately  0.1  bit/pixel  is  incurred  independently  of 
the  amount  of  compression,  to  keep  track  of  the  number 
of  components  truncated  per  line  and  the  scale  factors 
relating  the  dynamic  range  (both  real  and  imaginary) 
of  each  line  to  the  maximum  dynamic  range.  This 
overhead  is  included  in  the  results. 


Summary 

We  have  shown  that  the  Radon  transform  can  be  used 
to  advantage  in  bandwidth  compression  for  several 
reasons.  First,  a  line  passing  through  the  center  of  the 
two-dimensional  Fourier  transform  of  the  object  is  at¬ 
tainable  by  a  one-dimensional  operation  that  can  be 
carried  out  by  existing  fast  devices.  Second,  the  dy¬ 
namic  range  of  this  line  is  reduced  in  a  filtering  opera¬ 
tion  required  by  the  inverse  Radon  transform.  This 
reduction  enhances  compression  performance  (i.e., 
quantization  error  is  reduced).  Finally,  each  line  in 
Fourier  space  is  obtained  independently,  so  the  com¬ 
pression  can  be  adapted  to  the  amount  of  structure  in 
that  line.  The  technique  was  demonstrated  on  a  rec¬ 
tilinear  object,  and  a  bit  rate  of  1.6  bits/pixel  was 
achieved  with  good  fidelity. 

This  research  was  sponsored  by  the  U.S.  Air  Force 
Office  of  Scientific  Research  under  grant  AFOSR-82- 
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